{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Finding the axes of a hyper-ellipse\n",
    "\n",
    "Let us try finding the axes of the hyper ellipse described by the\n",
    "<br>equation $5x^2+6xy+5y^2=20$.\n",
    "<br>Note: The actual ellipse we use as example is 2D (to facilitate \n",
    "<br>visualization), but the code we develop will be general\n",
    "<br>and extensible to multi-dimensions.\n",
    "\n",
    "The ellipse equation can be written using matrices and vectors as\n",
    "<br>$\\vec{x}^{T}A\\vec{x} = 1$ where\n",
    "$A=\n",
    "\\begin{bmatrix} \n",
    "        5 & 3 \\\\\n",
    "        3 & 5 \\\\\n",
    "\\end{bmatrix} \n",
    "\\space \\space \\space\n",
    "\\vec{x} = \n",
    "\\begin{bmatrix} \n",
    "        x \\\\\n",
    "        y \\\\\n",
    "\\end{bmatrix}$.\n",
    "\n",
    "To find the axes of the hyper ellipse, we need to transform the\n",
    "<br>coordinate system so that the matrix in the middle becomes diagonal.\n",
    "<br>Here is how this can be done:\n",
    "<br>If we diagonalise $A$ into $S\\Sigma S^{-1}$, then the ellipse equation\n",
    "<br>becomes $\\vec{x}^{T}S \\Sigma S^{-1}\\vec{x} = 1$ where $\\Sigma$ is a\n",
    "<br>diagonal matrix.\n",
    "<br>Since $A$ is symmetric, its eigenvectors are orthogonal.\n",
    "<br>Hence, the matrix containing these eigenvectors as columns is orthogonal,\n",
    "<br> i.e., $S^{-1} = S^{T}$. In other words, $S$ is a rotation matrix.\n",
    "<br>\n",
    "<br>So the ellipse equation becomes $\\vec{x}^{T}S \\Sigma S^{T}\\vec{x} = 1$\n",
    "<br>or $\\left(\\vec{x}^{T}S\\right) \\Sigma \\left(S^{T}\\vec{x}\\right) = 1$\n",
    "<br>or $\\vec{y}^{T} \\Sigma \\vec{y} = 1$ where $\\vec{y} = S^{T}\\vec{x}$.\n",
    "<br>This is of the desired form since $\\Sigma$ is a diagonal matrix.\n",
    "<br>Remember, $S$ is a rotation matrix. Thus, rotating the coordinate system\n",
    "<br>by $S$ aligns the coordinate axes with the ellipse axes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 5 x^{2} + 6 x y + 5 y^{2} = 20$"
      ],
      "text/plain": [
       "Eq(5*x**2 + 6*x*y + 5*y**2, 20)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import torch\n",
    "from sympy import Symbol\n",
    "import sympy as sy\n",
    "\n",
    "x = Symbol('x')\n",
    "y = Symbol('y')\n",
    "a = Symbol('a')\n",
    "b = Symbol('b')\n",
    "ellipse_eq = sy.Eq(5*x**2 + 5*y**2 + 6*x*y, 20)\n",
    "ellipse_eq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sujaynarumanchi/.virtualenvs/book_test_pyt113/lib/python3.8/site-packages/sympy/plotting/plot.py:916: MatplotlibDeprecationWarning: \n",
      "The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.\n",
      "  self.ax.spines['left'].set_smart_bounds(True)\n",
      "/Users/sujaynarumanchi/.virtualenvs/book_test_pyt113/lib/python3.8/site-packages/sympy/plotting/plot.py:917: MatplotlibDeprecationWarning: \n",
      "The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.\n",
      "  self.ax.spines['bottom'].set_smart_bounds(False)\n",
      "/Users/sujaynarumanchi/.virtualenvs/book_test_pyt113/lib/python3.8/site-packages/sympy/plotting/plot.py:956: MatplotlibDeprecationWarning: \n",
      "The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.\n",
      "  self.ax.spines['left'].set_smart_bounds(False)\n",
      "/Users/sujaynarumanchi/.virtualenvs/book_test_pyt113/lib/python3.8/site-packages/sympy/plotting/plot.py:957: MatplotlibDeprecationWarning: \n",
      "The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.\n",
      "  self.ax.spines['bottom'].set_smart_bounds(False)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAGJCAYAAADv+MuDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deZzNdf/G8dfYZiwzYxs7NZbSRpIlZElZ0qJQiEQpa4pK2rRYSoo7IZKlEJUtWmVXobKEUESEQjFjyzbn98f7xxhnMOv5nOV6Ph7n4WwzrnvuMV3z+X6WMI/H40FERETkLFlcBxARERH/o4IgIiIiXlQQRERExIsKgoiIiHhRQRAREREvKggiIiLiRQVBREREvKggiIiIiBcVBBEREfGigiAiIiJeVBBEMki7du1YvHix6xgiIhlCBUEkg8TFxXHzzTdTrlw5BgwYwM6dO11HShGPx0N8fDw6lkVEzqaCIJJBZs6cyc6dO+ncuTNTp07l0ksvpXHjxnzyySecOHHCdbzzOnjwINHR0Rw8eNB1FBHxIyoIIhkoJiaGnj17smbNGpYvX07ZsmVp27YtxYoV4/HHH+e3335zHVFEJEVUEEQywe7du5k7dy5z584la9as3Hrrraxdu5Yrr7ySIUOGuI4nInJRYR5deBTJECdOnODTTz9l3LhxfP3111SoUIGHHnqI1q1bExUVBcCMGTPo0KED+/fvd5w2UXx8PNHR0cTFxZ3JKSKSzXUAkWBRtGhREhISaNWqFStWrODaa6/1ek+9evXImzevg3QiIqmjEQSRDPLBBx/QokULIiIiXEdJFY0giEhyVBBEQpwKgogkR5MURURExIsKgoiIiHhRQRAREREvKggiIiLiRQVBREREvKggiIiIiBcVBBEREfGigiAiIiJeVBBERETEiwqCiIiIeNFhTX5s1y44dOj8rz/9NPzxR8b/vVFRMGrUhd8TGQlFi2b83y0iIv5BBcFPvP46HDmS9LmpU2HDBjd5Lr/8wq9fdRW0aOH9fMmS0KFD5mQSERHf0WFNPvTzz3DypN2fOxfGjk18bfNmSEhI++eOioIyZdKX73y2b4d//knZe8PD4ZJLkj43dCgULpz4uEgRKFYs4/JJ+uiwJhFJjgpCJlq2DKZNS3w8fDgcPZq2z1W4MPTqdf7Xy5aFu+5K2+e+mK+/hjVrLvye115LeYmoXh1q1Up8nC8fPPNM2vNJ+qggiEhyVBAy0OrV8O+/0LatPT58GOLiUvax1apBjhx2v0EDaN8+6evZsiX9Ldzf/P134ujI2fr3h3Xrkr7v11+TvidLFhtVOK1zZ6hd2+5HREDVqhmfVxKpIIhIclQQ0mnpUpg0ye5Pnw579qTs42rVgtatEx+3awe5cmV8Pn+zaRPMn5/0uS5dzv/+PHmgTZvExw88YGVKMo4KgogkRwUhDebPh23b4PHH4fhx+O+/87/3sssS5waMHAl589r9HDkgZ85MjxoQzh1lGTECliyx+0eOwKJFia9FRNjXLkuWxGIG9nUuWzbzswYjFQQRSY4KQiq88w58+63NK7jYXIKrr4annoJKley+pM3BgzBzZuLjd96B777zfl/FinDNNXb/rbdsXoOkjAqCiCRHBeECDh2Czz+HhQth9Gg4der8761a1WbvDx9u/3EKC4OsWX0WNWScOgWnv2P/+gt69rT7a9fCxo12P0sW+/pXqGB7RYBNjCxVyvd5A4EKgogkRwXhHMeP23LD7t1tE6K5c8//3vBwqFcPHn4Yrr/e9gAQN9avT5z8ePJk4kTRkyetVFx/PZQoAdmzw/vv22tZs9rjUKeCICLJUUE4y6pV0LUrfP/9hd/XogWULg2vvuqbXJJ2CxcmFoLZs2HfvsTX6tZNLBK33w4xMb5O5x9UEEQkOSFfEE6dsnkFr7wCW7bA1q3e78mRw7YW7tIFbrzRbhERvs8q6bNihU2I3LkTnnjCnjt82CaZVqkChQrB+PH2fESEraAIBSoIIpKckC8IL74IL72U/GuRkYnL6u67z5epxFdmzYJ58+z+u+8mrkipVStxK+mWLa08BCsVBBFJTkgWhH37bPLa+vW22+G5Che23yCnT4frrvN9PnHjhx9szsLKlXY2xr//2iqKq6+2yw9jx9rSVH/esCotVBBEJDkhVxCWLoUePew/Aslp0sSWyZUu7dtc4n8+/zxxP4YhQ+DYMShfHu6/3567++6LH2oVCFQQRCQ5IVMQDh+GCRPsksLevUlfy5ULYmOhWTN47DGtoRdvv/xi81V++gkGD4bduyF3boiOhjFj7H5kpPdBVYFABUFEkhMSBeHAAVu2OHGi92vR0TBsWOJsdpGUmDvXJj2+/bbtxwBQrpzNW2jSBGrUcJsvNVQQRCQ5QV8QVq2Cpk3tyOJzVawIn36qDXQk7XbutL0ztm2DPn1sZOHgQSuer75ql6qqVHGd8sJUEEQkOUFdEIYPt2OE4+O9X+vSBQYMsB/kIhll1SpYvNjunz4Cu1MnK6MdOrjNdj4qCCKSnKAsCFu3Qu/eMGOG9xHEl14KgwbBnXcmHq8skhn+/df2XejY0S5z/fGHPf/kk7Y6pm5dO8bbNRUEEUlO0BWEZctsN8TkVilUr27XjCtX9n0uCW1bt9opoJMmwYIF9lzbtjZv4fnn3WZTQRCR5ARNQTh61JYnPv88nDiR9LWICHj0UXj5ZTs/QcSV48dtVOupp+DPP22EYflym9h4331WYosU8W0mFQQRSU7QFIS77kp6LPDZZsywiYoi/uaff+Cbb+yciMmTbQfHokVty+fwcDuZMrOpIIhIcnzw4yfzzZplqxHOddNNNqtc5UD8VYECcO+9tgQ3IcEOjcqe3c6BKFfOnl+1ynVKEQlFAT2CcPSonaMwbBgcOZL4fPbstu9Bnz5QsKC7fCJp9cUXdnz188/bMdWxsbY9eMWKtiFTWFjG/V0aQRCR5AR0QXjkERg92vv5Z56B/v19n0ckM3z5pV0+++wzm7fw3nsZu2RSBUFEkuMHi6xS78gR2zb53MsK2bLZ5YSHHnKTSyQzNGpkt59/hv374cEHrSxcdpmNlEVH21bPIiIZKSBHED78EFq3Tvpc7twwcKD9wBQJdtu22aW15ctt1cMNN9ioQlrOEdEIgogkJ+AKws8/2yZH27Ylfb5mTTupUSSU7N5t83Dmz7eSXLCg7SAaGWmrIVJCBUFEkhNQqxgOHbK9Ds4tB127wsKFLhKJuFW0KLzzjk1ofOABuPpqO4K6YUMrDSIiaRVQIwjffAO33JL0uTJlYN68wDxmVySjnTgBv/8OS5bAlCk2X6dTJysOlSolv/pBIwgikpyAGkG47Tbv50aMUDmQ0DZw4ECqVKlCZGQkxYsXonfvptx44ya++QZeeAH+/tuKtVb2iEhqBExBqF0bjh1L+tzAgfa8SChbtGgRXbt2ZdmyZcydO5cTJ07QoEEDDh8+TKNGdjjU+vV2yaFxY2jZ0k6c/Osv18lFxJ8FxCWGt9+G556zfetPa9AA3nwTrrrKXS4Rf7R3714KFSrEokWLqH1Og961y3Ye/ewzOHwYypeH3r3jiY3VJQYRScrvRxDWrYNx45KWgyJFoFUrlQOR5MT9/z+W/Pnze71WrBh06HCMsWPjqVs3ns2b47nuungAfvjBpzFFxM/5/QhC5842S/tsLVvaXggiklRCQgJ33HEHBw4cYOl51v2++OKLvPTSS17PFywYR//+UTz8cGanFJFA4NcFYdAg6N076XPVqsG330LWrG4yifizzp0788UXX7B06VJKlCiR7HuOHTvGsbMm9MTHx1OyZEmWLo3jrbei2LDB5vbUrm3HUGuXRpHQ5LdbLZ86BcePJ30uOhpee03lQCQ53bp1Y86cOSxevPi85QAgPDyc8PBwr+evuQamToUNG2DjRtuAadw4Oy49SxbIkSMz04uIv/HbOQjTptlJdmdr2xbq1HGTR8RfeTweunXrxowZM5g/fz6xsbHp+nxXXAF33QWrV9s8nzZt7BTJsWPtOREJDX55iWHjRmjRwiYonu3oUYiIcJNJxF916dKFyZMnM2vWLC6//PIzz0dHR5MzZ86LfvzFNko6csQu6y1YALNnw8svw4032rkPGs0TCV5+WRBWrLC5Bme77z4b7sye3U0mEX8Vltz2iMC4ceN44IEHLvrxqdlJ8bPP4Kuv4OOPbYvz555LS2IRCQR+VxD+/dcmR61fn/hcTAx8/71tqywiGSstWy2vWwft2kFsLOTPDwMG2EFRIhI8/G4OwsmTScsBwIQJKgci/uTqq+Gnn2wJ8sGD0KiRlYRz/+2KSODyu4KwaVPSx0WKQJ48brKIyIUVLAijR9vlv8OH4d57Yfp0FQWRYOB3BaFhw6SPX3jBJkSJiH+KjLQlkv37Q79+tkyyZk0YMsR1MhFJD78rCCISuJo2hWeftcsP06fb4VAPPQQ7d7pOJiKp5VcFYeDApJsj3XST/SYiIoGlTBlYssTmKMTH27LlHj1g2zbXyUQkpfyqIEycaDsonnb99VChgrs8IpI+l1xik4xnzYKoKNvobMYM16lEJCX8dqtlEQkOOXPa7ZVXoEoVWx45Z47NW3jsMdfpROR8/GYfhLvugpkzkz6eOlUbI4lktrTsg5AeR4/Cn3/CAw9AgQLQoYMVh+LFM/2vFpFU8JtLDOcezJQtm8qBSDDKmRPKlbPdGO+7z854aNECfvsNEhJcpxOR0/ymIJwtf34tbRQJdsWK2b4Jn34KLVva6ZH58tl5DyLinl/MQTh6FE6cSHxcpgx07+4uj4j41qOP2p+zZ0PHjrY08qmnIDzcbiLie34xgvDEEzB3rusUIuLa7bfbHgr79sGtt0Lr1vD7765TiYQmvygIIiKnFS9u+ydMnWqXG1u0gKVLYfdu18lEQovfFYSICJu4JCKhrXhxePddu/wwcyaULw/jx7tOJRI6/GIOwtkiI23HNRERsH0TAJo1g/bt4bvv7LJkTIxNahSRzOF3IwgiIsm54QZYuBCKFrVJzK1a2cFQIpI5nBeEv/+GvXtdpxCRQFCkiC2HnDIFKle2UYXp02HNGtfJRIKP84LwwQe2YYqISErly2fHSw8YAFu22BkPOl5aJGM5LwgiImnVtCk8+ST8+KONKtx5pxUGEUk/vysI/fq5TiAigaZsWTsAqkgRO8elc2dddhBJL78rCG3auE4gIoEoJgaGDYPFi21JZNOmtp/C2bu0ikjK+d0yRxGRtMqRw249ekCFCnD//Xb5YcwY18lEAo/fjSCIiGSEevVg61YbQWjUyJZJrl7tOpVI4FBBEJGglS0bTJgAH35oKx1atYJZs+DUKdfJRPyfCoKIBL18+eDVV20p5NNPQ8+esH+/61Qi/k0FQURCRqNGMGMG/PKLzU9o1gwOHHCdSsQ/qSCISEgpX96Ol58925ZFVqgAn3/uOpWI/9EqBhEJWUOH2uTFnj1h1So74yEqynUqEf+gEQQRCVnZs9veKx99BN98A23bwsCBsG+f62Qi7qkgiEjIq1DBSkK/frB+PdStCz//7DqViFsqCCIi2E6M11wDEydCt25w4422E6NIqFJBEBE5R6dOtmXz6NHQujVs3+46kYjvqSCIiCSjYkX47DObp3D33bZ984YNrlOJ+I7zglCqFBQv7jqFiIi3okVh1Cj44gtbEtmkCbz/PiQkuE4mkvmcF4R77rEz3EVE/FFEhM1P6NMHRoyAJ5+0ZZEiwc55QRARCRSNGsEff8CWLbby4ccfXScSyTx+VxC6dXOdQETk/CIibElkixa2VfOHH4LH4zqVSMbzu50U58xxnUBE5MJy5oTnn4f69aFzZ9tkqV49uOMO7cQowcPvRhBERAJFjRrwySe2ymHkSGjaVIc/SfBQQRARSYdy5Wx1w7ffQs2aULo0zJzpOpVI+vlFQXjsMahWzXUKEZH0eeUVeO89W+Xw7LMQH+86kUja+UVBKFcOChRwnUJEJP3uugvmzYOVK+3Sw9q1rhOJpI1fFISznTypbU1FJLDFxtqZDtWq2VbN06fD8eOuU4mkjt8VhP37oXlz1ylERNKnQAHo399OiOzVyzZY0nJICSR+UxDKlIHwcNcpREQy1p132sFPmzbZJMZff3WdSCRl/KYgvPWWlQQRkWBTsiR8+aXtlXDTTTBliutEIhfnNwVBRCTYvfQSjBsHr70G3bvDX3+5TiRyfn5ZELZtg3ffdZ1CRCRjZcsGt9xiIwhbt0KrVrrkIP7LrwrCsmW2z/nevTB3rus0IiKZ4/LLbQfGRo1su+YJE1wnEvHmV2cxREZCWJjrFCIimS8iAnr3hkqVbCnkmjXw5puuU4kk8qsRhLPt2GH/YEREglmDBvD337bK4ZprYMUK14lEjN8VhFGj7M9ly2DyZLdZRER8IWtWO0L63nttXsKkSXDqlOtUEur8riA0bOg6gYiI7+XODc89Z/MRXn8dHn5YqxzELb8rCGc7cgSOHnWdQkTEd2rVsgmM//5r5zps3Og6kYQqvysIOXNCixZ2/+23Yfx4p3FERHyubFmYMQNatoSqVbXKQdzwu4IQGQlPPOE6hYiIez162MmQL78MDz0E+/a5TiShxO8KgoikzuLFi7n99tspVqwYYWFhzJw503UkyUBVqtjKhiNHbN+EpUtdJ5JQ4ZcF4eqrbZIOwCuvwJIlbvOI+LPDhw9TsWJFhg8f7jqKZJICBWD0aGjbFjp1gpEj4fBh16kk2PnVRkmn5coFpUrZ/d274dAht3lE/Fnjxo1p3Lix6xiSyfLksUsOFSvan99/D8OGQXS062QSrPxyBAEgf34oUcLub9ig1QwiGeXYsWPEx8cnuUngqFvXtqIPC7PLD6tXu04kwcpvC8LNN0O3bna/Vy/44w+3eUSCxcCBA4mOjj5zK1mypOtIkkqFCtnKhl69rDCMGOE6kQQjvy0IIpI5+vTpQ1xc3Jnbjh07XEeSNHr4YVi82MrCvffaCZEiGcWvC0KjRlCnjt1/9lk4ccJtHpFgEB4eTlRUVJKbBKawMKhQAaZNg/BwaN4cli93nUqChV8XhIoV4Yor7P706dqbXEQkOSVKwNixdlm2ZUsYNAj++891Kgl0frmK4WwPPmizddesUUEQSc6hQ4fYvHnzmcdbt25l9erV5M+fn1KnlwNJ0MuWDdq3h/Ll4YEHbCTho4/sICiRtPDrEQSA66+HokXtfpkybrOI+KMff/yRSpUqUalSJQB69uxJpUqVeOGFFxwnExduuMGOji5UCK66ynUaCWR+XxDA9kUIC4OEBG0OInKuunXr4vF4vG7jdZBJSBs6FPbvhz59tExc0iYgCsK0aXDZZbB3ry1/FBGRCwsPt0uza9ZAgwawapXrRBJoAqIgAMTEuE4gIhJYihSBiRNtr4T777dfto4fd51KAkXAFIQlS6wR79oFc+a4TiMiEhjy57czbV59FZ55xjZX0qVaSYmAKQgAsbGwfTtMmuQ6iYhIYGnSBL74AjZvhvr1Yds214nE3wVUQViwwP5cuxa++cZtFhGRQFO6tJWEO+6AW291nUb8XUAVhGzZoFw5WL8e5s1znUZEJDA9/jjkyGHbM2/f7jqN+KuAKggFC9qe4wDz58OyZW7ziIgEopw5bRQ2d24bTVi40HUi8UcBVRAAIiPh8sthxQpbviMiIqlXsKCdAtmjB3TpYtszHzzoOpX4k4ArCFdfDX372v1Jk+CXX9zmEREJVBERtj3zmDEwZYpt0XzokOtU4i8CriCAbb182WW29PHPP12nEREJbDVq2M/TLVtsEqMIBGhBqFvXmi7Aa6/B7t0u04iIBL7cuaF/f3j4YdueWSQgCwLYiWWxsTZZUdfNRETSr0kTG5VduRKqVLE/JXQFbEG46y5o3NgOcerSBU6ccJ1IRCTw5c4NH38MDRtCq1Ywa5brROJKNtcB0mP4cPjtN5g7Fzwe12lERIJDVBT06weFC8PkyVCxIlx6qetU4msBO4JwWlSUjSK0bes6iYhIcLnzTlvp0KIF/PADJCS4TiS+FPAF4ZNPoGxZ+PJLGxYTEZGMUaqUbU73yCO2NfPQoa4TiS8FfEEAKF4c4uP1zSsikhkeesi2tx83TqO1oSQoCsKCBZA9O+zYoQk1IiKZoUIFWLrU5nxptDY0BEVBANvcY8cO+Ogj10lERIJTnjy26+JLL8GTT8K+fa4TSWYKmoJQtSp07AirV8PXX7tOIyISfLJmhdtug+nTbZv7Vq203X0wC5qCEBkJ7drZN+vkybBnj+tEIiLB6bLL7OdsiRIwezYcPuw6kWSGoCkIYDsrtmljs24//dR1GhGR4BUdbVveT5xoo7dHj7pOJBktqApCsWL2DVu6NDz2GCxa5DqRiEjwqlPHJi7mzAlXXQXLl7tOJBkpqAoCQP360KGDDXlp+2URkcwVHQ3vvQdPPQWdOrlOIxkp6AoC2PWxUqVsY48tW1ynEREJfm3aQL58cOONsHat6zSSEYKyILRoAQ0a2AhC9+7aHlREJLPlyQPTpkG1atC8uS09l8AW0Ic1Xci778LWrfZNqoOcREQyX758MHiwzQebPBmuvhpKlnSdStIqKEcQTuvXD4oUgfbtXScREQkdTZvCf//ZaO66dfolLVAFdUGoXt32EJ8zB2bOdJ1GRCQ0lC5t2zHffz/Urg2HDrlOJGkR5vEEf7crWtSGub76yobARCRRfHw80dHRxMXFERUV5TqOBJmbboLwcHj/fYiJcZ1GUiOoRxBOW7HCzjLv0cN1EhGR0DJ/PlxxBVx7rTawCzQhURCyZbMjoefMgZEjXacREQktgwbBa6/B8OGwd6/rNJJSIVEQihaFDz+E/fthyhRYtcp1IhGR0JEtm11qiI6GO++0UV3xfyFREAAqV4ZXX4XFi1UQRER8rVgx23GxWTNbWbZwoetEcjEhUxBy5bJjSm+5BT76CDZvdp1IRCS0REZCr15Qrx689Rb89ZfrRHIhIVMQwA4Tad3aRhBmz7Z1uiIi4lv9+tlR0dWrw5EjrtPI+YRUQQA77bFJE+jZE3bvdp1GRCT05M1rIwi5c8Mzz7hOI+cTcgUB4O237UCRChW0w5eIiCvLltnl3nr1XCeR5IRkQciVC557DiIi4PHHXacREQlNkZEwbhz8/DMMG+Y6jZwrJAsC2GmP99wDEyfCggWu04iIhKaCBeG772wp+p13agK5PwnZggC2acehQ9C3rx0NLSIivhUWBpdfbqvLoqPh88/h1CnXqQRCvCAAbNoES5ZAx46uk4iIhK4SJaBVKxgxwpZCHjjgOpGEfEEoVAhef90Ocnr/fddpRERCV+PGMGMGfPst/PKLJpG7FvIFIWdOu+511VV2DWzDBteJRERC1xVXQIcOthz9nXdcpwltIV8QAMqVg8mTbSvQzz7TBkoiIi517myTx/v10341Lqkg/L9ChWyXxQkTYPly12lERELbtddCy5ZQtapddhDfU0E4S/36cPvt0LWrlj6KiLj26qswcCCMGqVzG1xQQTjHgAFQoIAdTSoiIu5kz267LEZHw91365hoX1NBSMaUKTYfoXZt10lEREJb8eJ2THTNmjBtmpY/+pIKQjKKFoWXX4bVq+3cBhERcSdPHtv5dskSO3Bvzx7XiUKDCsJ5dOgAV18NY8bYelwREXGnShXbryYuznZbTEhwnSj4qSCcR1gYfPkl3HwzDB7sOo2IiERGQu/eduvZ03Wa4KeCcAFRUXD//XbS2MCBrtOIiEijRvDHH7ax3eLFrtMENxWEi6hQwZY9Tp4MCxdq608REdciImyeWO/eutSQmVQQUqB9ezvxsX17WLrUdRoREbnnHhvlrVXLJpRLxlNBSKHmzW0+Qp8+sHev6zQiSQ0fPpxLL72UiIgIqlWrxgotGJcgly+fLUmvUQMmToTjx10nCj4qCKnw7rs2efG991wnEUk0depUevbsSd++fVm5ciUVK1akYcOG7NFaMAly+fJBixY2F2HCBNdpgk+Yx6Or6qmxejU0aGCjCZMnu04jAtWqVaNKlSq8/f+bdiQkJFCyZEm6d+/O008/fdGPj4+PJzo6mri4OKKiojI7rkiG+9//YNw4eO012zI/WzbXiYJDir6MHo+HgwcPZnaWgFC6NDz0kB1DOn68bf8p4srx48f58ccf6dGjB/Hx8Weer127NosXL6ZLly5eH3Ps2DGOHTt25vHpf9tnf7xIIGnfHvLmhQcfhNmzoUwZ14n8T2RkJGFhYan6mBSNIJz+DUNEREQCT1pGCFNUEDSC4O2PP+KpUKEkV121g9dfj6JmTdeJglN8fDwlS5Zkx44dGv5Oxu7duylfvjxz586latWqZ55//vnn+fbbb5k/f77Xx5w7grB7926qVq3KL7/8QvHixX2SO9To+9g3Tn+dixTZwaZN+jqfLS0jCCm6xBAWFqZv6nNccon9eeedUcyZE0W1apA/v9tMwSwqKkrfg8mIiIgga9asHDp0KMnX58CBAxQvXjxVX7PIyEh9jTOZvo99o3DhKDp1itI8sXTSKoZ0uvVW2LABvvjCdRIJRTly5KBy5crMmzfvzHMJCQnMmzePG264wWEyEXcmT4ZFi2y3RUk7FYR0uuYaaNXKdvUaMwbOGrkV8YmePXvy7rvvMmHCBDZs2EDnzp05fPgw7du3dx1NxIkiRWDYMFvZsG2b6zSBS4tB0ig8PJy+ffsSHh5Ox44QGwsPP2wjCsWKuU4XPM7+Okvy7r33Xvbu3csLL7zAX3/9xbXXXsuXX35J4cKFU/Txp7+2+hpnHn0f+8bpr3Pu3OFUrw6ffgrNmllZqFHDdbrAo30QMlCvXjaktWuX6yQiKad9ECRYzZsHnTpBkybw3HNQsKDrRIFFlxgy0BtvQJYs8OSTrpOIiEj9+jBtGqxcaafySuqoIGSwLVtg3Tpo3Nh1EhERqVDB5om1bm2XGjRmnnIqCBksPBzeegvWroW334aTJ10nEhEJbZ07Q8WK8OijmkieGioImaBcOfjyS1tq8/rrrtOIiMhXX9lk8ssvd50kcKggZJKrr4ZHHrGS8OmncOKE60QiIqHt/fdhzx6bUC4Xp4KQTtu2bePBBx8kNkobraUAAB5tSURBVDaWnDlzUqZMGfr27cvx48dp1w4GDLBvxp9+cp008PXv358aNWqQK1cu8ubN6zpO0Bg9ejQAMTExVKtWjRUrVjhOFFwWL17M7bffTrFixQgLC2PmzJmuIwWdgQMHUqVKFSIjIylUqBBNmzZl06ZNXu+rVQtatoSZM2HuXAdBA4wKQjpt3LiRhIQERo0axfr16xkyZAjvvPMOzzzzDAC33w6NGkHXrrBmjeOwAe748eO0aNGCzp07u44SNKZOnXrme3XJkiVUrFiRhg0bsmfPHsfJgsfhw4epWLEiw4cPdx0laC1atIiuXbuybNky5s6dy4kTJ2jQoAGHDx/2eu8779gvbaNGwYEDDsIGEo9kuEGDBnliY2OTPDdihMeTP7+jQEFm3LhxnujoaNcxgkLVqlU9HTt29ACeuLg4z6lTpzzFihXzDBw40HW0oAR4ZsyY4TpG0NuzZ48H8CxatCjZ13//3eNp0sTj6dvXt7kCjUYQMkFcXBz5zzm5qXNnm5fQsKG2/hT/cPz4cX766Sfq1q175rksWbJw88038/3337sLJpJOcXFxAF4/h0+LjYWHHoJPPoHhwyGZgQZBlxgy3ObNmxk2bBiPPPKI12vz5sEVV1hZEHFt3759nDp1ikKFCiV5vnDhwvz111+OUomkT0JCAo899hg1a9bk6quvPu/7mjaF5s2hWzf44QcfBgwgKgjn8fTTTxMWFnbB28aNG5N8zM6dO2nUqBEtWrSgY8eOXp8zWzbo0QP++gs6dtSWzJC2r7OIyPl07dqVdevWMWXKlIu+t0cPuPtum7j4558+CBdgdFjTefTq1YsHHnjggu8pXbr0mfu7du2iXr161KhR48ys8OTExsKMGbYd88SJ8NRTGZU4MKX26ywZp2DBgmTNmtVrQuLff/9NkSJFHKUSSbtu3boxZ84cFi9eTIkSJS76/nz5bDn6mjXwwgswdqwPQgYQFYTziImJISYmJkXv3blzJ/Xq1aNy5cqMGzeOLFkuPDBz6aXQti08/bTt6vXooxAdnQGhA1Bqvs6SsXLkyEHlypVZtGjRmecSEhKYN28e3bp1c5hMJHU8Hg/du3dnxowZLFy4kNjY2BR/bIMGcOedNhehQwdbCilGBSGddu7cSd26dbnkkksYPHgwe/fuPfPahX4Lu+MOO7O8Sxe7/jVxIuggvQvbvn07//77L9u3b+fUqVOsXr0agLJly5InTx7H6QJTz549uf/++wHYtGkTY8aM4fDhw7Rv395xsuBx6NAhNm/efObx1q1bWb16Nfnz56dUqVIOkwWPrl27MnnyZGbNmkVkZOSZOTTR0dHkzJnzoh//xhs2YbFtW9i6NbPTBhDXyygC3bhx4zxAsreUOHbM46lc2eP54INMDhoE2rVrl+zXecGCBa6jBbRBgwZ5AE/27Nk9VatW9Sxbtsx1pKCyYMGCZL9v27Vr5zpa0Djfz+Bx48al+HMcOeLxgMfTsGHm5Qw0YR6PzrZybe5cO22sXTtrsiK+FB8fT3R0NHFxcURpGEtClMdjB+y98gq89JJWm4FWMfiFW26BP/6AKVNg/nzXaUREQk9YGNx7rx0P/eGH8PPPrhO5p4LgJ8LD4fnnoX17GDnSdRoRkdBTqJD9HF6xwn5hO37cdSK3VBD8RLZs0KkTjBgBU6fC/8+/ExERH6pTx34ODxxoJSGUqSD4merVoUYNuP9+2y8h1BusiIiv5cwJuXND//6hfalBBcHPFChgR0T37w+9e4f2N6eIiAutWkHPnvDrr3b6Y3y860RuqCD4qdtvt4OdHnoIVq1ynUZEJLR06gS33mpzwv7/7KeQo4Lgx4YNs42U6tWza2IiIuIbxYrZL2glS8LQoa7TuKF9EALA6tV26thXX0GZMq7TSLDRPggi5/fggzBpEvz0E1x1les0vqURhABw5ZW2NXODBjBzpus0IiKh4733ICbGzmsINSoIASBHDnjtNRgyxE4cu8BhkSIiksGWLbMzGu6913US31JBCBDZs9sowkMPwfjxmrgoIuIrRYpAQgIsXQoffOA6zcXt3buXIkWKMGDAgDPPfffdd+TIkYN58+al+POoIASYRx6B1q3h7rtt4uKJE64TiYgEtyxZYOVK2LXLTt4963BOvxQTE8PYsWN58cUX+fHHHzl48CBt27alW7du1K9fP8WfR5MUA9TChXa40+rVkC+f6zQSyDRJUeTi4uNtf5pBg6BfP3j2WdeJLq5r16588803XH/99axdu5YffviB8PDwFH+8RhACVN26dj2sTBn4+mvXaUREgltUlB2sB7B2rR2w5+8GDx7MyZMn+fjjj5k0aVKqygGoIAS0QYNsl6/77rNmq8sNIiKZp0IFOwZ66lQ7eTchwXWiC9uyZQu7du0iISGBbdu2pfrjdYkhCGzfbqdAXnNN6G7oIWmnSwwiKbd4MXTvbtvgL10KNWu6TpS848ePU7VqVa699louv/xyhg4dytq1aylUqFCKP4dGEIJAqVLQq5dtpPT++3DsmOtEIiLBqXZtuOceW35+4ACcOuU6UfKeffZZ4uLieOutt+jduzeXXXYZHTp0SNXn0AhCEFm4EB57DK6/3o4qjYlxnUgCgUYQRFLvvvtg8mT4808oXtx1mqQWLlzILbfcwoIFC6hVqxYA27Zto2LFirz66qt07tw5RZ9HBSHIbN0KTzxhE2gWLoQ8eVwnEn+ngiCSej162Bywa6+F5ctdp8kcusQQZGJjYdo0m7D41luu04iIBKf//c9Ga3/80YpCMFJBCFKffAIzZtix0du3u04jIhJ8YmMhLAzefNN1ksyhghCkypWz2balStna3c8/d51IRCS4TJwIhQvDvn0wbpzrNBlPBSGI5cxpzfbll6F3b5gwwXUiEZHgsngx7N8PH37oOknGU0EIcuHhtuNi+/YwdiysW+c6kYhI8ChZ0v7cuNE2UAomKggholMnuO02O9N87Fg4edJ1IhGRwJc9O3z/PezYYfvQbN3qOlHGUUEIEblywZNPwsiR8Pzz8OijrhOJiAS+sDDIn9/uf/45TJniNk9GUkEIMQ0awJYt8Ndf/rtFqIhIIClbFj791O6vX2+bJwUDFYQQFBEBo0fbNbM33vD/A0dERPxZliw23wtg0iRYssRtnoyighCiChaEH36AmTNtbsKmTa4TiYgErgYNEjdMWrIE9u51mycjqCCEsNKl4aOPbK+Eli1h0SL/PXhERMTfZc9ut5EjYfNm12nSTwUhxBUtaq33qaegbVvbM+HECdepREQCT4cOdrIu2E62R464zZNeKggCQKtW8NlnMH588EywERHxtchIm+f1+usQH+86TfqoIMgZ11wDXbtCpUq2haiIiKTOM89A06Z2f/Rot1nSSwVBknjqKZu42K+flYVgmGgjIuJLRYrYqoYBA1wnSR8VBPFSty58843tldC0qa12EBGRlBkyxEZiT56Et992nSbtVBAkWSVK2OlkTZvCAw/Y/f/+c51KRCRwnDqlgiBBKirKtmceMgRefRW6dNEKBxGRlOjWDfLkgT17YMwY12nSRgVBLqpBA/j2W9sn4bvvXKcREfF/990HuXPbUdDTprlOkzYqCJIiBQvaxMXbb4dBg1ynERHxf/36uU6QPioIkmKtWsHatXYoSePGdpaDiIgkr00b+3PuXHj8cbdZ0kIFQVLlkktsh7DYWGjWDBYscJ0otPXv358aNWqQK1cu8ubN6zqOiCTj1KnAnOStgiCpFhNjM3Nffhk6doTnnoODB12nCk3Hjx+nRYsWdO7c2XUUETlHeLiNuoJtaR9ol2ezuQ4ggSlLFhtBKF0aHn4YVqywkYXcuV0nCy0vvfQSAOPHj3cbRES8hIVBtrP+K+vxuMuSFhpBkHSpVMk2UqpcGUqWhDlzXCeSizl27Bjx8fFJbiKSOaKj4ZZb7P4PP8DWrW7zpIYKgmSIgQNtCK1rV3jxRTh82HUiOZ+BAwcSHR195layZEnXkUSCVtGiiRMUp00LrJ1pVRAkw9xzj01a/O47uOsuWL/edaLA9PTTTxMWFnbB28Z0LCHp06cPcXFxZ247duzIwPQiEizCPJ5Auyoi/m7vXjvq9KuvbB1wo0aQPbvrVIFj7969/PPPPxd8T+nSpcmRI8eZx+PHj+exxx7jwIEDqf774uPjiY6OJi4ujqioqFR/vIhc2P79dnDT4MFQtizMmgVXXuk61cVpkqJkuJgYm61brRp07w7Nm9vjLBqvSpGYmBhiYmJcxxCRDJIvn03oBti8GQ4dcpsnpfQjWzJNs2Ywfz6sXAn16sHvv7tOFHy2b9/O6tWr2b59O6dOnWL16tWsXr2aQ4HyE0gkBG3ZYic9+jsVBMlUpUtbSaheHerUgY8/dp0ouLzwwgtUqlSJvn37cujQISpVqkSlSpX48ccfXUcTkbNUrQq1atn91q3h33/d5kkJzUEQnzh5Er7+Gvr0gZtvht69oVAh16kENAdBxFd6907cLOnvv/3/Z6BGEMQnsmWDW2+FKVPgl1+sQf/2m+tUIiJyPioI4lNXXAGffGKXG265BSZPDrzdxURE0qJsWZvEDXDddW6zpIQKgvhc7tzw/PMwYgQ8+qgNu4mIBLuOHaFuXbsfCHMQVBDEmVtvhR07YM0a26p55UrXiURE5DQVBHEqZ06YPh1uu82WRX70kS45iEjwCguzP48dg1693Ga5GK1iEL+xeDF06wY33GBHSRcu7DpRaNAqBhHfqlrVzmS48kr/3pJeIwjiN2rXtgmMf/4Jd99tm4mIiIgbKgjiVy67DD77zArCtdfCpEmuE4mIZI5//4XPP3ed4vxUEMQv9eqVuLFSly6QhjOIRET82l9/wfDhrlOcnwqC+K0bboBVq+x0yAYN4PvvXScSEUm/1q0hIsJ1iotTQRC/VqAAvPce3HOPrSF+7z04csR1KhGRtHvsMYiMdJ3i4lQQxO9FRcETT8D//me3bt0gPt51KhGR4KaCIAGjfn346is4etROh/Tn5UEiIimxcCG89JLrFMlTQZCAUrQofPghdO0KNWvC6NGuE4mIpN2RI7Bnj+sUyVNBkIDUpQvMmwcff2x7mz/+uP/+IxMRCUQqCBKQwsLs/Ia5c+3Qpw0bbGawNlcSEckYKggS8K68EmbOtJGEOnXsCGkREX/288+uE1xcNtcBRDJCRAQ89xxcf72NJKxaBa+/7jqViEjyAuHYE40gSFBp1Ah27YJ16+C66+DHH10nEhG5sM2bYdMm1ym8qSBI0ImIsGOj77gD2rSB5s1h/nxISHCdTETE29dfw+zZrlN4U0GQoBQZCS++CBs3QsuW8OCD8PzztoeCiIhcnAqCBL3mzWHWLFsWeeedsH+/60QiIv5PBUFCQoUKsGwZVKoE5ctbYRARkfNTQZCQMmCAHa/aty80aQLDhsHhw65TiYj4HxUECSlZs9olh9WrbU7CpEk2P2HXLtfJRET8iwqChKzq1WHGDMiTBxo0sEsQIiJiVBAkpBUtCmPGQK9etofCM8/AhAmuU4mIuKeCIAK0bw9r1lhJGDoUHnlElx1EJLSpIIj8v0sugdq1Yfp02LcP7rnH5iqIiIQiFQSRc8TGwsSJ0KyZTWgcPRoWLHCdSkTEt1QQRJKRMyc8/jiMGmUrH+6+25ZGioiEChUEkQuoX9+WQW7ZAkuXQs2agXFMq4hIeqkgiKRA/vwwbZoVhHvusfMdli0Dj8d1MhEJRKdOuU5wcSoIIimUNy8MGmQHQN12G7RoAf37w/HjrpOJSKApU8Z1gotTQRBJgzZtbJOlqVPh3nshPh6OHHGdSkQk46ggiKTR9dfD2rVw6aXQoQPccAN8+63rVCISaK68EipWdJ3CmwqCSDoNGQLjx0PbttCxI7z3nm2ydOKE62QiEgjq1oVbbnGdwpsKgkgGyJMHnngC3nkHNmyAG2+E3r3h2DHXyURE0kYFQSQD1a4NgwfDV1/ZcshbboFff9W2zSKSaPt2OHnSdYqLU0EQyQRly8I330DjxvDGG7afwldfuU4lIv6geXPYv991iovL5jqASDDr0wf++w8++sguQfz0kx0IdcUVtlujiIi/0giCSCaLiID774exYyFHDrjvPujcWfMTRATKl4e77nKdInkqCCI+UqWKjSIsXWqbK117LSxeDL/84jqZiPjKL7/YvimnlS4NN9/sLs+FqCCI+FiBAjB5MvTsCZ99ZrsyTp/uOpWI+MKLL8KmTa5TpIwKgogjHTvaVs0jR8L778M119jjgwddJxMR0SRFEaeyZYOGDe22fDk8/LAdAlWqFAwf7jqdiGSkKVPs3/lpxYpBr17u8lyMRhBE/ES1arByJcycCXv2wGWXwfffu04lIhll1SrbA+G0vHnhppvc5bkYFQQRP5I1q90++MAuQbRta6dGzp4dGMfDikjyTp703n49PNxNlpQK83h0or2IP/vsM5vQeOutNk+hXTsrERklPj6e6Oho4uLiiIqKyrhPLCJnjBwJXbokPs6Z0/9PgNUcBBE/16SJLYX688/ElQ+jRtlvH5GRrtOJSLDSJQaRAHDFFXauw9q1NoGxc2c7Ae6771wnE5GL2bYNFixI+twDD7hIkjq6xCASgA4etFGEGTMgVy7o3t3OfciePfWfS5cYRDLXF1/YJcKz/f03FCrkJk9KaQRBJABFRtqujN9+a9c1H33UjpceNsx1MhE52759trwxEGkOgkiAu+su27b50CFo0wYWLbLd2mJioHBh1+lEQldCgl0GfP/9pM9PmgT587vJlBq6xCASROLj4e23bb31rl3wv//B9ddf7GN0iUEkM/z3X/Knti5fDlWr+j5PamkEQSSIREXBM8/YaMLo0TaSsGMHvPwylCwJ113nOqFI6Fi1yvu5YcOgYkXfZ0kLzUEQCVDbtm3jwQcfJDY2lpw5c1KmTBn69u3L8ePHyZPHlkTOmQN9+9os6rp14a23XKcWCR3nHuN8223QoIH/b5B0mkYQRALUxo0bSUhIYNSoUZQtW5Z169bRsWNHDh8+zODBg8+87+677c8GDaw0TJtmZ9C3aGEjCtn0U0AkwzVubCsVzlarlm2hHig0B0EkiLz++uuMHDmS33///bzv+e03WLwYxo61SxIjRsRTurTmIIhklFGj4PnnYe/exOfq1LE5QYFyeQF0iUEkqMTFxZH/ItOjS5U6RosW8cyaFc+VV8bz6KPxgF2OWLnSFylFgtf69TBxYtJyEBMDLVsGVjkAFQSRoLF582aGDRvGI488csH3DRw4kOjoaGJionnzzWjmzCkJ2BLJ5s3ho4/sYCiNLYqkTkKCLWlcujTp81WqQKdObjKlhy4xiPiZp59+mtdee+2C79mwYQPly5c/83jnzp3UqVOHunXrMmbMmAt+7LFjxzh27NiZx/Hx8ZQsWZK4uDhWrYri3Xdt7Xbr1jZXoU2b9P3vEQkVo0Z5F4HLL4eNG93kSS8VBBE/s3fvXv75558Lvqd06dLkyJEDgF27dlG3bl2qV6/O+PHjyZIldQODye2DsHo17NxpWzjXqQNDh9o2zrlype1/k0iwW7kS7r/fLjGcli0bzJ1rK4gCkQqCSADbuXMn9erVo3LlykycOJGsaTgH+kIbJSUk2FbOf/4J//5rPwBvuimwZmKLZLYDB+CFF7y3Ou/QAd57z02mjKCCIBKgdu7cSd26dbnkkkuYMGFCknJQpEiRFH+elOykuH8//PyzzU/47jtb/fD223bYjLZzllDm8dicg9q1kz5foAD88ov/H8h0IVoBLRKg5s6dy+bNm9m8eTMlSpRI8lpG9/58+exSQ506MHmyDafecANERMDTT8M110DDhhn6V4oEhOHD7VLcuWbPDuxyABpBEAl5aT2L4eefYft223wpe3Zo1w6aNrUtnZPbf14k2GzfbmednL2kEaBfP+jRA/LkcZMro2iZo4ikSYUKtnXsr7/a0q4jR6BZM+jWDV55xXU6kcw1bRrUrOldDm66Ce64I/DLAWgEQSTkZeRpjrt2wcGD0LmzrYKoUAH69IFixSAV0yJE/NqQITYp8dChpM/fcw+MGGHzD4KBRhBEJMMUK2brvr/+Gjp2hLx5oXJlO/fhiSfgjz9cJxRJn+XLYcoU73JQrhw8+GDwlAPQCIJIyMvIEYRz/fcf/PMPfPghfPqpLZc8vdtcvnw2uVEkEBw9aueXPPssxMUlfa1FCxs5KFjQTbbMooIgEuIysyCca9o0KwmPPWaTGZs0sSNxGzTI1L9WJF0SEmxH0Q8/9H4te3Y4ftz3mXxBBUEkxPmyIICtG4+Ph1WrYNAg2LbN5iv06wdlytjkrnPXlIu48v33tu34tm3er1WrZqUhNtbnsXxCBUEkxPm6IJzrp59sQ5kXXrAfwlFRNgv85ptt6aSIK++8Y9+X565UADtz4aWXAn+vgwtRQRAJca4LwmmnTtmchfbtbTXEt99Clix2OaJaNbjzTggPdxZPQsjWrfDUU/DJJ96vXXopvP463H23fX8GMxUEkRDnLwXhbLt3w7Jldn/IEPjhB1tfXqYMnD7oMnt2OwxHJKMkJNguoR072oFl57r+ejux8brrfJ/NhSDvPyISiIoWtcmLd90FixfbDPLSpW3nuly57Namjc0qv8jBlyIp4vHYyECVKsmXg969raiGSjkAjSCIhDx/HEE4n4MHbWThiy9sqWR8vP3AjoyESy6BV1+18yGCYRc78Z24OHjkEZg61fu1qCgbNWjZ0ve5XFNBEAlxgVQQzjVmjJ0JcewYjB5tz9WpY9eHb7vNRh1Ezmf5cru99RZs2eL9+s0320TFMmV8n80fqCCIhLhALginnTxpw78//ABvvmnP5cuXOJrw7rv2XFQU5M/vLqf4h//+s9GCZ5+1Jbbnyp3bVin07g0xMb7P5y9UEERCXDAUhOS8/z5s3AiDB8OJE/ZczZpw66228125cm7ziRsDBtgqhTFjkn+9VCkrmc2a+TaXP1JBEAlxwVoQTlu/HjZvtt8WwZZQRkbaaEJkZOKliYIFdaBUsPrzTzhwwI5gnj8/+fdcdpkta8ydW5emTlNBEAlxwV4QzvXll7Y5E8CLL9rlCbC9FurXt/s9egT3BjihYtUq+4/+nDk2V+V8OnSwUYPoaN9lCwQqCCIhLtQKwtm2bbPRheees8dbt8KePXYqZY4cULYsvPKKvRYbC4ULO4sqKbRqlY0SdetmJy7u25f8+woWtM2Qmje3kaOcOX2bMxCoIIiEuFAuCOeaPx/WrrX7TzyROLoAULcuVKhg9994Q5s0+ZspU+zchPHjbfnr+eTNayNHV1yhQ8IuRgVBJMSpICTv9P77330Hw4bZ/ZUrYf9+KFAAwsLsuaZN4d577X6NGraJk2S+Y8dgyRJYuhSGD7c9Mo4dS/695ctD8eIwdKhtwlWggG+zBioVBJEQp4KQcl9+acviDhywEYZzNWuWeB37hhvgoYd8my8UvPmmTTw9fDj5jY3OVr26/X9w4402CVFSRwVBJMSpIKSex2Nr6QFmzICPPkp8bdYs+zNLFpvHcNrQoYlzGEqWhMqVfZM1kP33n5UysK/rlCk2SnCh/2rVr29LFYcPh6xZk/5/IKmjgiAS4lQQMtZHH9l/wJYsgREjEp8/+yftJZfYqomzffhh0sdhYYmXMYKZx5P0a/Pee/DNN3b/6FGYPfvCH3/6a9Sxox3odeuttnxV0k8FQSTEqSD4xuOPJ85r2LIl8bTK83n4YRsaP61IEdv6NxisWAG//mr34+Oha9fUf46WLW1uQd++GZtNEqkgiIQ4FQTf+/NPWLMm8XH37kmX4/33X+Luj6flz2/X1M9VpIgdiX0h2bNn3jK+I0eSrvY411NPwY4dSZ9bt85O5kyJs0cDHnwwsSQ1bmyXcSTzqCCIhDgVBP/z6afw+eeJj+Pi7Pp7WtWqBffdl/5cyRk58sKbEKVWrVpw1VV2PyoKBg3KuM8tqaOCIBLiVBD839GjNix/thdfhE2bvN+7f3/iBEp/Ex7ufVjW2dtdg52RUayYb3NJ8lQQREKcCkJwmTjRdhO8kJEjrXRktJgYaNv2/K9XqADt2mX83yuZQwVBJMSpIISeNWsuPG8grXLmhCuvzPjPK25os1ARkRBTsaLrBBIINAdUREREvKggiIiIiBcVBBEREfGigiAiIiJeVBBERETEiwqCiIiIeFFBEBERES8qCCIiIuJFBUFERES8qCCIiIiIFxUEERER8aKCICIiIl5UEERERMSLCoKIiIh4UUEQERERLyoIIiIi4iXM4/F4XIcQEXc8Hg8HDx4kMjKSsLAw13FExE+oIIiIiIgXXWIQERERLyoIIiIi4kUFQURERLyoIIiIiIgXFQQRERHxooIgIiIiXlQQRERExMv/AQYJFfUI/TByAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note that the ellipse major axis is forming an\n",
      "angle of 45 degrees with X axis.\n",
      "Rotating coordinate system by this angle\n",
      "will align ellipse axes with coordinate axes.\n"
     ]
    }
   ],
   "source": [
    "# Let us plot this ellipse\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy as sy\n",
    "\n",
    "def sym_eq_plot(eq):\n",
    "    plot = sy.plot_implicit(eq)\n",
    "    return plot \n",
    "plot = sym_eq_plot(ellipse_eq)\n",
    "\n",
    "print(\"Note that the ellipse major axis is forming an\\n\"\n",
    "      \"angle of 45 degrees with X axis.\\n\"\n",
    "      \"Rotating coordinate system by this angle\\n\"\n",
    "      \"will align ellipse axes with coordinate axes.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigen values are:  tensor([8.+0.j, 2.+0.j])\n",
      "\n",
      "Eigen vectors are columns of S matrix\n",
      "tensor([[ 0.7071+0.j, -0.7071+0.j],\n",
      "        [ 0.7071+0.j,  0.7071+0.j]])\n",
      "\n",
      "Rotation angle theta = 45.00 degrees\n"
     ]
    }
   ],
   "source": [
    "A = torch.tensor([[5, 3], [3, 5]]).float()\n",
    "\n",
    "# Obtain eigen values and vectors of the ellipse\n",
    "# coeeficients matrix\n",
    "l, S = torch.linalg.eig(A)\n",
    "\n",
    "print(\"Eigen values are:  {}\\n\".format(l))\n",
    "print(\"Eigen vectors are columns of S matrix\\n{}\".format(S))\n",
    "\n",
    "l = l.real\n",
    "S = S.real\n",
    "\n",
    "# Assert that eigen vectors are orthogonal\n",
    "assert torch.dot(S[:, 0], S[:, 1]) == 0.0\n",
    "\n",
    "# Find the angle between the principal axis and the X-axis.\n",
    "import math\n",
    "\n",
    "# Vector corresponding to X-axis\n",
    "x_axis_vec = torch.zeros((A.shape[0]))\n",
    "x_axis_vec[0] = 1\n",
    "\n",
    "# First principal eigen vector\n",
    "first_eigen_vec = S[:, 0]\n",
    "\n",
    "# Dot product between the two vectors (equals cosine\n",
    "# of the angle between the directions of the two vectors)\n",
    "dot_prod = torch.dot(x_axis_vec, first_eigen_vec)\n",
    "\n",
    "# The angle between the two vectors is the cosine inverse\n",
    "# of the dot-product, in radians\n",
    "theta = math.acos(dot_prod)\n",
    "\n",
    "# Convert to degrees from radian\n",
    "theta = math.degrees(theta)\n",
    "print(\"\\nRotation angle theta = {:.2f} degrees\".format(theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGdCAYAAADuR1K7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3de3RU5aH38d+EkATFJOWWMZAItlQQKbRgQjjnPVgzy6B0CRWXmIWAlEqpgNZQCigXL+1JFS+ggBzPqaVUKTTU0kopHhq8UBm5BFRu4WirXJ0ExEwQJYnJ8/6RZnQkCQlmz+WZ72etvdA9z848z16RfJ3Ze+IyxhgBAABYIi7cEwAAAGhLxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAq8SHewLhUFdXp+PHj+uSSy6Ry+UK93QAAEALGGN0+vRppaenKy6u6ddnYjJujh8/royMjHBPAwAAXIAjR46oR48eTT4ek3FzySWXSKo/OcnJyWGeDQAAaInKykplZGQEfo43JSbjpuGtqOTkZOIGAIAoc75LSrigGAAAWIW4AQAAViFuAACAVYgbAABgFeIGAABYhbgBAABWIW4AAIBViBsAAGAV4gYAAFiFuAEAAFYhbgAAgFWIGwAAYBXiBgAAWIW4AQAAViFuAACAVYgbAABgFeIGAABYhbgBAABWIW4AAIBViBsAAGAV4gYAAFiFuAEAAFYhbgAAgFWIGwAAYBXiBgAAWIW4AQAAViFuAACAVYgbAABgFeIGAABYhbgBAABWIW4AAIBViBsAAGAV4gYAAFiFuAEAAFYhbgAAgFWIGwAAYBXiBgAAWIW4AQAAViFuAACAVYgbAABgFeIGAABYhbgBAABWIW4AAIBVQhI3S5cuVc+ePZWUlKTs7Gxt37692fFFRUXq06ePkpKS1L9/f23YsKHJsVOmTJHL5dKiRYvaetoAACAKOR43a9asUUFBgRYsWKBdu3ZpwIABysvLU3l5eaPjt27dqvz8fE2aNEm7d+/WqFGjNGrUKO3du/ecsX/84x/1xhtvKD093ellAACAKOF43Dz++OO64447NHHiRF155ZVavny5LrroIj377LONjl+8eLGGDx+umTNnqm/fvnrooYf0ne98R0uWLAkad+zYMU2fPl3PP/+82rdv7/QyAABAlHA0bqqrq1VSUiKPx/P5E8bFyePxyOv1NnqM1+sNGi9JeXl5QePr6uo0btw4zZw5U/369TvvPKqqqlRZWRm0AQAAOzkaNydPnlRtba3S0tKC9qelpcnn8zV6jM/nO+/4hx9+WPHx8brrrrtaNI/CwkKlpKQEtoyMjFauBAAARIuou1uqpKREixcv1ooVK+RyuVp0zJw5c+T3+wPbkSNHHJ4lAAAIF0fjpkuXLmrXrp3KysqC9peVlcntdjd6jNvtbnb8li1bVF5erszMTMXHxys+Pl6HDh3SjBkz1LNnz0a/ZmJiopKTk4M2AABgJ0fjJiEhQYMGDVJxcXFgX11dnYqLi5WTk9PoMTk5OUHjJWnTpk2B8ePGjdPbb7+tN998M7Clp6dr5syZeumll5xbDAAAiArxTj9BQUGBJkyYoMGDBysrK0uLFi3SmTNnNHHiREnS+PHj1b17dxUWFkqS7r77bg0bNkyPPfaYRowYodWrV2vnzp165plnJEmdO3dW586dg56jffv2crvduuKKK5xeDgAAiHCOx82YMWN04sQJzZ8/Xz6fTwMHDtTGjRsDFw0fPnxYcXGfv4A0dOhQrVq1SnPnztW9996r3r17a926dbrqqqucnioAALCAyxhjwj2JUKusrFRKSor8fj/X3wAAECVa+vM76u6WAgAAaA5xAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqIYmbpUuXqmfPnkpKSlJ2dra2b9/e7PiioiL16dNHSUlJ6t+/vzZs2BB4rKamRrNmzVL//v118cUXKz09XePHj9fx48edXgYAAIgCjsfNmjVrVFBQoAULFmjXrl0aMGCA8vLyVF5e3uj4rVu3Kj8/X5MmTdLu3bs1atQojRo1Snv37pUkffLJJ9q1a5fmzZunXbt26YUXXtDBgwd14403Or0UAIhKdXXSX/8q1daGeyZAaLiMMcbJJ8jOztbVV1+tJUuWSJLq6uqUkZGh6dOna/bs2eeMHzNmjM6cOaP169cH9g0ZMkQDBw7U8uXLG32OHTt2KCsrS4cOHVJmZuZ551RZWamUlBT5/X4lJydf4MoAIPIVF0uzZkkjR0rz5oV7NsBX09Kf346+clNdXa2SkhJ5PJ7PnzAuTh6PR16vt9FjvF5v0HhJysvLa3K8JPn9frlcLqWmpjb6eFVVlSorK4M2ALDZ7t1SXp7k8UhHj0r33BPuGQGh42jcnDx5UrW1tUpLSwvan5aWJp/P1+gxPp+vVePPnj2rWbNmKT8/v8mKKywsVEpKSmDLyMi4gNUAQOT75z+lsWOl73xH+t//rd+3YIHUsWN45wWEUlTfLVVTU6NbbrlFxhg9/fTTTY6bM2eO/H5/YDty5EgIZwkAzisvl+66S+rTR1q16vP93/iG9MMfhm9eQDjEO/nFu3Tponbt2qmsrCxof1lZmdxud6PHuN3uFo1vCJtDhw5p8+bNzb73lpiYqMTExAtcBQBEro8/lh5/XFq4sP6fv+w//1Nq3z708wLCydFXbhISEjRo0CAVFxcH9tXV1am4uFg5OTmNHpOTkxM0XpI2bdoUNL4hbN555x397W9/U+fOnZ1ZAABEqJoaadky6etfr3/bqbGwufpq6eabQz83INwcfeVGkgoKCjRhwgQNHjxYWVlZWrRokc6cOaOJEydKksaPH6/u3bursLBQknT33Xdr2LBheuyxxzRixAitXr1aO3fu1DPPPCOpPmxuvvlm7dq1S+vXr1dtbW3gepxOnTopISHB6SUBQNjU1UlFRdLcudK77zY/9uGHJZcrNPMCIonjcTNmzBidOHFC8+fPl8/n08CBA7Vx48bARcOHDx9WXNznLyANHTpUq1at0ty5c3Xvvfeqd+/eWrduna666ipJ0rFjx/TnP/9ZkjRw4MCg53r55Zd1zTXXOL0kAAiLhtu6S0rOPzYvT/rud52fExCJHP+cm0jE59wAiDZLlkgFBfVvR7XE7t3Sl/7/D4h6EfE5NwCAtjFtWv0dUb/5jTR8ePNjx44lbBDbiBsAiBKpqVJ+vnTJJU2Pad9eeuih0M0JiETEDQBEiZqa+ldlioqaHnPnnVKvXqGbExCJiBsAiAJNhc2/7rWQVP+Kzn33hXZeQCQibgAgwjUVNlOn1t9B1XDD6c9+JnXtGvr5AZHG8VvBAQAXrrmweeqp+s+xueYaad8+fjkm0IC4AYAI1ZKwkaRbbqn/cL+LLw79HIFIRNwAQARqadhI0pgxhA3wRcQNAESY1oSNVH+LOIDPcUExAESQ1oYNgHMRNwAQIQgboG0QNwAQAQgboO0QNwAQZoQN0LaIGwAII8IGaHvEDQCECWEDOIO4AYAwIGwA5xA3ABBihA3gLOIGAEKIsAGcR9wAQIgQNkBoEDcAEAKEDRA6xA0AOIywAUKLuAEABxE2QOgRNwDgEMIGCA/iBgAcQNgA4UPcAEAbI2yA8CJuAKANETZA+BE3ANBGCBsgMhA3ANAGCBsgchA3APAVETZAZCFuAOArIGyAyEPcAMAFImyAyETcAMAFIGyAyEXcAEArETZAZCNuAKAVCBsg8hE3ANBChA0QHYgbAGgBwgaIHsQNAJwHYQNEF+IGAJpB2ADRh7gBgCYQNkB0Im4AoBGEDRC9iBsA+BLCBohuxA0AfAFhA0Q/4gYA/oWwAexA3ACACBvAJsQNgJhH2AB2IW4AxDTCBrAPcQMgZhE2gJ2IGwAxibAB7EXcAIg5hA1gN+IGQEwhbAD7ETcAYgZhA8QG4gZATCBsgNhB3ACwHmEDxBbiBoDVCBsg9hA3AKxF2ACxibgBYCXCBohdIYmbpUuXqmfPnkpKSlJ2dra2b9/e7PiioiL16dNHSUlJ6t+/vzZs2BD0uDFG8+fP16WXXqoOHTrI4/HonXfecXIJAKIIYQPENsfjZs2aNSooKNCCBQu0a9cuDRgwQHl5eSovL290/NatW5Wfn69JkyZp9+7dGjVqlEaNGqW9e/cGxjzyyCN68skntXz5cm3btk0XX3yx8vLydPbsWaeXAyDCETYAXMYY4+QTZGdn6+qrr9aSJUskSXV1dcrIyND06dM1e/bsc8aPGTNGZ86c0fr16wP7hgwZooEDB2r58uUyxig9PV0zZszQT3/6U0mS3+9XWlqaVqxYoVtvvfW8c6qsrFRKSor8fr+Sk5PbaKUAwo2wAezW0p/fjr5yU11drZKSEnk8ns+fMC5OHo9HXq+30WO8Xm/QeEnKy8sLjH/vvffk8/mCxqSkpCg7O7vJr1lVVaXKysqgDYB9ysulbduC9xE2QOxxNG5Onjyp2tpapaWlBe1PS0uTz+dr9Bifz9fs+IY/W/M1CwsLlZKSEtgyMjIuaD0AIlv37tLLv/qnMuOOSpKm9inWU4vrCBsgxsTE3VJz5syR3+8PbEeOHAn3lAA44ehRXT7u3/RK3f/TAt2vp0o9ct3xQ6muLtwzAxBCjsZNly5d1K5dO5WVlQXtLysrk9vtbvQYt9vd7PiGP1vzNRMTE5WcnBy0AbDQpZdK112nXnpf9+sBuSTp17+WfkjgALHE0bhJSEjQoEGDVFxcHNhXV1en4uJi5eTkNHpMTk5O0HhJ2rRpU2B8r1695Ha7g8ZUVlZq27ZtTX5NADGiXTvp2Wel8eOD9xM4QEyJd/oJCgoKNGHCBA0ePFhZWVlatGiRzpw5o4kTJ0qSxo8fr+7du6uwsFCSdPfdd2vYsGF67LHHNGLECK1evVo7d+7UM888I0lyuVz6yU9+op///Ofq3bu3evXqpXnz5ik9PV2jRo1yejkAIl1D4EjSypWf7//1r+v//J//keJi4h15IGY5HjdjxozRiRMnNH/+fPl8Pg0cOFAbN24MXBB8+PBhxX3hL5qhQ4dq1apVmjt3ru6991717t1b69at01VXXRUY87Of/UxnzpzR5MmTVVFRoX//93/Xxo0blZSU5PRyAEQDAgeIaY5/zk0k4nNugBhRWyv94AfBgSNJEycSOEAUiojPuQGAsOIaHCAmETcA7EbgADGHuAFgPwIHiCnEDYDYQOAAMYO4ARA7CBwgJhA3AGILgQNYj7gBEHsIHMBqxA2A2ETgANYibgDELgIHsBJxAyC2ETiAdYgbACBwAKsQNwAgETiARYgbAGhA4ABWIG4A4IsIHCDqETcA8GUEDhDViBsAaAyBA0Qt4gYAmkLgAFGJuAGA5hA4QNQhbgDgfAgcIKoQNwDQEgQOEDWIGwBoKQIHiArEDQC0BoEDRDziBgBai8ABIhpxAwAXgsABIhZxAwAXisABIhJxAwBfBYEDRBziBgC+KgIHiCjEDQC0BQIHiBjEDQC0FQIHiAjEDQC0JQIHCDviBgDaGoEDhBVxAwBOIHCAsCFuAMApBA4QFsQNADiJwAFCjrgBAKcROEBIETcAEAoEDhAyxA0AhAqBA4QEcQMAoUTgAI4jbgAg1AgcwFHEDQCEA4EDOIa4AYBwIXAARxA3ABBOBA7Q5ogbAAg3AgdoU8QNAEQCAgdoM8QNAEQKAgdoE8QNAEQSAgf4yogbAIg0rQ2cDz+UqqpCNz8gwhE3ABCJWhM4q1dLy5eHdn5ABCNuACBStTRwfv976aGHJL8/9HMEIhBxAwCR7HyBc+yYtGVL/VtTjz4anjkCEYa4AYBI11zgXHONZEz9vz/+uPTBByGfHhBpiBsAiAZNBc67737+z598Ij34YGjnBUQg4gYAokVTgfNF//3f0v/9X+jmBEQg4gYAooXPJ/3Xf0nvv9/0mNpaae7ckE0JiETx4Z4AAKAFliyR7r67ZR/iV1Qkbd8uZWU5Py8gAvHKDQBEg2nTpNdfl/7jP1o2ftaszy80BmKMY3Fz6tQpjR07VsnJyUpNTdWkSZP08ccfN3vM2bNnNXXqVHXu3FkdO3bU6NGjVVZWFnj8rbfeUn5+vjIyMtShQwf17dtXixcvdmoJABBZhgyRXnlF+stfpP79mx/7yivSSy+FYlZAxHEsbsaOHat9+/Zp06ZNWr9+vV577TVNnjy52WPuuecevfjiiyoqKtKrr76q48eP66abbgo8XlJSom7duum5557Tvn37dN9992nOnDlasmSJU8sAgMjickk33CDt3i395jdSZmbTY2fN4ndRISa5jGn71y0PHDigK6+8Ujt27NDgwYMlSRs3btQNN9ygo0ePKj09/Zxj/H6/unbtqlWrVunmm2+WJJWWlqpv377yer0aMmRIo881depUHThwQJs3b27x/CorK5WSkiK/36/k5OQLWCEARIizZ6Vly6Rf/EI6dercx3/7W+m220I/L8ABLf357cgrN16vV6mpqYGwkSSPx6O4uDht27at0WNKSkpUU1Mjj8cT2NenTx9lZmbK6/U2+Vx+v1+dOnVqu8kDQDRJSpIKCqR//lO6916pQ4fgx+fN45dqIuY4Ejc+n0/dunUL2hcfH69OnTrJ5/M1eUxCQoJSU1OD9qelpTV5zNatW7VmzZrzvt1VVVWlysrKoA0ArJKSUv/qzbvvSpMn138mjlR/2zi/VBMxplVxM3v2bLlcrma30tJSp+YaZO/evRo5cqQWLFig6667rtmxhYWFSklJCWwZGRkhmSMAhFx6ev1n4ezbJ40eXb+PX6qJGNOquJkxY4YOHDjQ7Hb55ZfL7XarvLw86NjPPvtMp06dktvtbvRru91uVVdXq6KiImh/WVnZOcfs379fubm5mjx5sua24MOq5syZI7/fH9iOHDnSmmUDQPS54gpp7VrJ65X69eOXaiKmtOpD/Lp27aquXbued1xOTo4qKipUUlKiQYMGSZI2b96suro6ZWdnN3rMoEGD1L59exUXF2v0v/5v4+DBgzp8+LBycnIC4/bt26drr71WEyZM0C9+8YsWzTsxMVGJiYktGgsAVmm4ffy11+o/vbjh7SrAYo7cLSVJ119/vcrKyrR8+XLV1NRo4sSJGjx4sFatWiVJOnbsmHJzc7Vy5Upl/etTNH/84x9rw4YNWrFihZKTkzV9+nRJ9dfWSPVvRV177bXKy8vTwoULA8/Vrl27FkVXA+6WAgAg+rT057djv37h+eef17Rp05Sbm6u4uDiNHj1aTz75ZODxmpoaHTx4UJ988klg3xNPPBEYW1VVpby8PC1btizw+Nq1a3XixAk999xzeu655wL7L7vsMr3f3O9aAQAAMcOxV24iGa/cAAAQfcL6OTcAAADhQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArOJY3Jw6dUpjx45VcnKyUlNTNWnSJH388cfNHnP27FlNnTpVnTt3VseOHTV69GiVlZU1OvbDDz9Ujx495HK5VFFR4cQSAABAFHIsbsaOHat9+/Zp06ZNWr9+vV577TVNnjy52WPuuecevfjiiyoqKtKrr76q48eP66abbmp07KRJk/Stb33LiakDAIAo5jLGmLb+ogcOHNCVV16pHTt2aPDgwZKkjRs36oYbbtDRo0eVnp5+zjF+v19du3bVqlWrdPPNN0uSSktL1bdvX3m9Xg0ZMiQw9umnn9aaNWs0f/585ebm6qOPPlJqamqL51dZWamUlBT5/X4lJyd/xdUCAIBQaOnPb0deufF6vUpNTQ2EjSR5PB7FxcVp27ZtjR5TUlKimpoaeTyewL4+ffooMzNTXq83sG///v168MEHtXLlSsXFtWz6VVVVqqysDNoAAICdHIkbn8+nbt26Be2Lj49Xp06d5PP5mjwmISHhnFdg0tLSAsdUVVUpPz9fCxcuVGZmZovnU1hYqJSUlMCWkZHRyhUBAIBo0aq4mT17tlwuV7NbaWmpU3PVnDlz1LdvX912222tPs7v9we2I0eOODRDAAAQbvGtGTxjxgzdfvvtzY65/PLL5Xa7VV5eHrT/s88+06lTp+R2uxs9zu12q7q6WhUVFUGv3pSVlQWO2bx5s/bs2aO1a9dKkhouF+rSpYvuu+8+PfDAA41+7cTERCUmJrZojQAAILq1Km66du2qrl27nndcTk6OKioqVFJSokGDBkmqD5O6ujplZ2c3esygQYPUvn17FRcXa/To0ZKkgwcP6vDhw8rJyZEk/eEPf9Cnn34aOGbHjh36wQ9+oC1btujrX/96a5YCAAAs1aq4aam+fftq+PDhuuOOO7R8+XLV1NRo2rRpuvXWWwN3Sh07dky5ublauXKlsrKylJKSokmTJqmgoECdOnVScnKypk+frpycnMCdUl8OmJMnTwaerzV3SwEAAHs5EjeS9Pzzz2vatGnKzc1VXFycRo8erSeffDLweE1NjQ4ePKhPPvkksO+JJ54IjK2qqlJeXp6WLVvm1BQBAICFHPmcm0jH59wAABB9wvo5NwAAAOFC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAqxA3AADAKsQNAACwCnEDAACsQtwAAACrEDcAAMAqxA0AALAKcQMAAKxC3AAAAKsQNwAAwCrEDQAAsApxAwAArELcAAAAq8SHewLhYIyRJFVWVoZ5JgAAoKUafm43/BxvSkzGzenTpyVJGRkZYZ4JAABordOnTyslJaXJx13mfPljobq6Oh0/flyXXHKJXC5XuKcTdpWVlcrIyNCRI0eUnJwc7ulYi/McGpzn0OA8hwbnOZgxRqdPn1Z6erri4pq+siYmX7mJi4tTjx49wj2NiJOcnMx/PCHAeQ4NznNocJ5Dg/P8ueZesWnABcUAAMAqxA0AALBKu/vvv//+cE8C4deuXTtdc801io+PyXcqQ4bzHBqc59DgPIcG57n1YvKCYgAAYC/elgIAAFYhbgAAgFWIGwAAYBXiBgAAWIW4iQGnTp3S2LFjlZycrNTUVE2aNEkff/xxs8ecPXtWU6dOVefOndWxY0eNHj1aZWVljY798MMP1aNHD7lcLlVUVDixhKjgxHl+6623lJ+fr4yMDHXo0EF9+/bV4sWLnV5KxFm6dKl69uyppKQkZWdna/v27c2OLyoqUp8+fZSUlKT+/ftrw4YNQY8bYzR//nxdeuml6tChgzwej9555x0nlxAV2vI819TUaNasWerfv78uvvhipaena/z48Tp+/LjTy4h4bf39/EVTpkyRy+XSokWL2nra0cXAesOHDzcDBgwwb7zxhtmyZYv5xje+YfLz85s9ZsqUKSYjI8MUFxebnTt3miFDhpihQ4c2OnbkyJHm+uuvN5LMRx995MQSooIT5/lXv/qVueuuu8wrr7xi/vGPf5jf/va3pkOHDuapp55yejkRY/Xq1SYhIcE8++yzZt++feaOO+4wqamppqysrNHxr7/+umnXrp155JFHzP79+83cuXNN+/btzZ49ewJjfvnLX5qUlBSzbt0689Zbb5kbb7zR9OrVy3z66aehWlbEaevzXFFRYTwej1mzZo0pLS01Xq/XZGVlmUGDBoVyWRHHie/nBi+88IIZMGCASU9PN0888YTTS4loxI3l9u/fbySZHTt2BPb99a9/NS6Xyxw7dqzRYyoqKkz79u1NUVFRYN+BAweMJOP1eoPGLlu2zAwbNswUFxfHdNw4fZ6/6M477zTf/e53227yES4rK8tMnTo18O+1tbUmPT3dFBYWNjr+lltuMSNGjAjal52dbX70ox8ZY4ypq6szbrfbLFy4MPB4RUWFSUxMNL/73e8cWEF0aOvz3Jjt27cbSebQoUNtM+ko5NR5Pnr0qOnevbvZu3evueyyy2I+bnhbynJer1epqakaPHhwYJ/H41FcXJy2bdvW6DElJSWqqamRx+MJ7OvTp48yMzPl9XoD+/bv368HH3xQK1eubPYXmMUCJ8/zl/n9fnXq1KntJh/BqqurVVJSEnSO4uLi5PF4mjxHXq83aLwk5eXlBca/99578vl8QWNSUlKUnZ3d7Hm3mRPnuTF+v18ul0upqaltM/Eo49R5rqur07hx4zRz5kz169fPmclHmdj+iRQDfD6funXrFrQvPj5enesjx9IAAAQ8SURBVDp1ks/na/KYhISEc/4CSktLCxxTVVWl/Px8LVy4UJmZmc5MPoo4dZ6/bOvWrVqzZo0mT57cNhOPcCdPnlRtba3S0tKC9jd3jnw+X7PjG/5szde0nRPn+cvOnj2rWbNmKT8/P2Z/AaRT5/nhhx9WfHy87rrrrrafdJQibqLU7Nmz5XK5mt1KS0sde/45c+aob9++uu222xx7jkgQ7vP8RXv37tXIkSO1YMECXXfddSF5TqAt1NTU6JZbbpExRk8//XS4p2OVkpISLV68WCtWrJDL5Qr3dCIGv6giSs2YMUO33357s2Muv/xyud1ulZeXB+3/7LPPdOrUKbnd7kaPc7vdqq6uVkVFRdCrCmVlZYFjNm/erD179mjt2rWS6u8+kaQuXbrovvvu0wMPPHChS4so4T7PDfbv36/c3FxNnjxZc+fOvbDFRKEuXbqoXbt259yp19g5auB2u5sd3/BnWVmZLr300qAxAwcObMvpRw0nznODhrA5dOiQNm/eHLOv2kjOnOctW7aovLw86BX02tpazZgxQ4sWLdL777/ftouIFuG+6AfOarjQdefOnYF9L730UosudF27dm1gX2lpadCFru+++67Zs2dPYHv22WeNJLN169Ymr/q3mVPn2Rhj9u7da7p162Zmzpzp3AIiWFZWlpk2bVrg32tra0337t2bvQDze9/7XtC+nJyccy4ofvTRRwOP+/1+Lihu4/NsjDHV1dVm1KhRpl+/fqa8vNyZiUeZtj7PJ0+eDPq7eM+ePSY9Pd3MmjXLlJaWOreQCEfcxIDhw4ebb3/722bbtm3m73//u+ndu3fQLcpHjx41V1xxhdm2bVtg35QpU0xmZqbZvHmz2blzp8nJyTE5OTlNPsfLL78c03dLGePMed6zZ4/p2rWrue2228wHH3wQ2GLpB8Xq1atNYmKiWbFihdm/f7+ZPHmySU1NNT6fzxhjzLhx48zs2bMD419//XUTHx9vHn30UXPgwAGzYMGCRm8FT01NNX/605/M22+/bUaOHMmt4G18nqurq82NN95oevToYd58882g79+qqqqwrDESOPH9/GXcLUXcxIQPP/zQ5Ofnm44dO5rk5GQzceJEc/r06cDj7733npFkXn755cC+Tz/91Nx5553ma1/7mrnooovM97//ffPBBx80+RzEjTPnecGCBUbSOdtll10WwpWF31NPPWUyMzNNQkKCycrKMm+88UbgsWHDhpkJEyYEjf/9739vvvnNb5qEhATTr18/85e//CXo8bq6OjNv3jyTlpZmEhMTTW5urjl48GAolhLR2vI8N3y/N7Z98b+BWNTW389fRtwY4zLmXxdLAAAAWIC7pQAAgFWIGwAAYBXiBgAAWIW4AQAAViFuAACAVYgbAABgFeIGAABYhbgBAABWIW4AAIBViBsAAGAV4gYAAFiFuAEAAFb5/4Tt5cO/rcsoAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the eigen vectors\n",
    "plt.quiver([0, 0], [0, 0], S[:,0], S[:,1],\n",
    "           color=['r','b','g'], scale=5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAAGiCAYAAAB+sGhNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzdd1gU58IF8LN0EQELVUFRo2hEwIICRjCiYMfeUWOPJrHEQm4siUmIRqNXYyyxd8XYO6JYELGBUSxRo6EIWNkVIqAw3x9+coMCAjI7u7Pn9zzzJDv7zu5hxT1OVwiCIICIiEjG9KQOQEREJDaWHRERyR7LjoiIZI9lR0REsseyIyIi2WPZERGR7LHsiIhI9lh2REQkeyw7IiKSPZYdERHJnqhlFxISgqZNm6JChQqwtrZGYGAgbt68+c7lQkND4ezsDBMTE7i4uODAgQNixiQiIpkTtexOnDiBMWPG4OzZswgLC8OLFy/Qtm1bZGRkFLrMmTNn0LdvXwwdOhQxMTEIDAxEYGAgrl69KmZUIiKSMYU6LwT98OFDWFtb48SJE2jZsmWBY3r37o2MjAzs27cvb17z5s3h5uaGpUuXqisqERHJiIE630ypVAIAKlWqVOiYqKgoTJgwId88f39/7Nq1q8DxWVlZyMrKynucm5uLJ0+eoHLlylAoFGWQmoiI1EkQBDx79gz29vbQ0yubDZBqK7vc3FyMGzcO3t7eaNCgQaHjUlJSYGNjk2+ejY0NUlJSChwfEhKCb775pkyzEhGR9BISElCtWrUyeS21ld2YMWNw9epVnD59ukxfNzg4ON+aoFKphKOjIxISEmBubl6m70VEROJTqVRwcHBAhQoVyuw11VJ2Y8eOxb59+3Dy5Ml3trStrS1SU1PzzUtNTYWtrW2B442NjWFsbPzWfHNzc5YdEZEWK8tdUaIejSkIAsaOHYudO3fi2LFjcHJyeucynp6eCA8PzzcvLCwMnp6eYsUkIiKZE3XNbsyYMdi0aRN2796NChUq5O13s7CwQLly5QAAQUFBqFq1KkJCQgAAX3zxBXx8fDBv3jx06NABW7ZswYULF7B8+XIxoxIRkYyJuma3ZMkSKJVK+Pr6ws7OLm/aunVr3pj4+HgkJyfnPfby8sKmTZuwfPlyuLq6Yvv27di1a1eRB7UQEREVRa3n2amDSqWChYUFlEol99kREWkhMb7HeW1MIiKSPZYdERHJHsuOiIhkj2VHRESyx7IjIiLZY9kREZHsseyIiEj2WHZERCR7LDsiIpI9lh0REckey46IiGSPZUdERLLHsiMiItlj2RERkeyx7IiISPZYdkREJHssOyIikj2WHRERyR7LjoiIZI9lR0REsseyIyIi2WPZERGR7LHsiIhI9lh2REQkeyw7IiKSPZYdERHJHsuOiIhkj2VHRESyx7IjIiLZY9kREZHsseyIiEj2WHZERCR7LDsiIpI9lh0REckey46IiGRP1LI7efIkOnXqBHt7eygUCuzatavI8REREVAoFG9NKSkpYsYkIiKZE7XsMjIy4OrqisWLF5douZs3byI5OTlvsra2FikhERHpAgMxX7xdu3Zo165diZeztraGpaWlCImIiEgXaeQ+Ozc3N9jZ2aFNmzaIjIwscmxWVhZUKlW+iYiI6N80quzs7OywdOlS/P777/j999/h4OAAX19fXLp0qdBlQkJCYGFhkTc5ODioMTEREWkDhSAIglreSKHAzp07ERgYWKLlfHx84OjoiPXr1xf4fFZWFrKysvIeq1QqODg4QKlUwtzc/L0yExGR+qlUKlhYWJTp97io++zKgoeHB06fPl3o88bGxjA2NlZjIiIi0jYatRmzILGxsbCzs5M6BhERaTFR1+zS09Nx+/btvMd3795FbGwsKlWqBEdHRwQHByMpKQnr1q0DACxYsABOTk748MMPkZmZiRUrVuDYsWM4cuSImDGJiEjmRC27CxcuoFWrVnmPJ0yYAAAYNGgQ1qxZg+TkZMTHx+c9n52djYkTJyIpKQmmpqZo2LAhjh49mu81iIiISkptB6ioixg7NomISH3E+B7X+H12RERE74tlR0REsseyIyIi2WPZERGR7LHsiIhI9lh2REQkeyw7IiKSPZYdERHJHsuOiIhkj2VHRESyx7IjIiLZY9kREZHsseyIiEj2WHZERCR7LDsiIpI9lh0REckey46IiGSPZUdERLLHsiMiItlj2RERkeyx7IiISPZYdkREJHssOyIikj2WHRERyR7LjoiIZI9lR0REsseyIyIi2TOQOgDJjyAIuPHoBv56+hfSs9PzpowXGUjPToeZkRmqmVdD1QpVX/3XvCrMjMykjk1EMsayo/f2utwi7kUg4u8IRNyLwIOMB/nGmBqaorxheZQ3Ko9nWc/w+PnjfM9bGFvA0cIRHlU90MKxBbwdvFG7Um0oFAp1/ihEJFMsOyo1QRCw7899CA4PRtzDOBjoGcCjqgeGug+Fbw1fNLBugApGFWBqaAp9Pf18yz5/8Rz3n91HoioRSc+SkKhKxJ0ndxCVGIVVMasgQIB1eWt4OXjB28EbLRxboIl9Exjo8VeWiEpOIQiCIHWIsqRSqWBhYQGlUglzc3Op48hWVEIUphydglPxp9CqRitM9p6Mjxw/Qnmj8u/92mmZaTibeBaR8ZE4nXAa0YnReP7yOaqYVkFX567o9WEv+NbwZfERyZQY3+MsOyqR5y+eY/je4dh4ZSNcrF0w2282AmoHiLq58UXOC1xKvoSdN3Yi9Foo/nr6F4uPSMZYdsXAshNPWmYaOm/ujAv3L2Bx+8UIcg16a/Ok2ARBwKXkSwi9FvpW8Q1xG4Lm1ZpzPx+RlmPZFQPLThwZ2Rlova41bj25hX1998HTwVPqSPmKb2vcVtxLuwd3W3eM9RiLvg36opxhOakjElEpsOyKgWVX9nJyc9Bpcyecij+FiEERaGzfWOpIb8kVcnHkzhH8cu4XHLh1ABXLVcRQ96EY3WQ0nCo6SR2PiEpAjO9xUU8qP3nyJDp16gR7e3soFArs2rXrnctERESgUaNGMDY2Ru3atbFmzRoxI1IxbLyyEQdvH0Roz1CNLDoA0FPoIaB2APb124dbn93CELch+O3Sb6i1sBY6b+6MI3eOIFfIlTomEUlE1LLLyMiAq6srFi9eXKzxd+/eRYcOHdCqVSvExsZi3LhxGDZsGA4fPixmTCpC5stMTDs+Dd3rdUdA7QCp4xRLrUq1MLftXCRNSMLyTssRr4yH/wZ/1FtcD6tjVuNFzgupIxKRmqltM6ZCocDOnTsRGBhY6JgpU6Zg//79uHr1at68Pn36IC0tDYcOHSrW+3AzZtmad2YephydgrhP41C3Sl2p45SKIAiITIjEvKh52HVjF2pY1sBU76kY7DYYxgbGUscjojdo3WbMkoqKioKfn1++ef7+/oiKiip0maysLKhUqnwTlY20zDR8f+p7DG80XGuLDnj1D60Wji2ws/dOXB51GR5VPTB6/2jUWlgLv5z7BVkvs6SOSEQi06iyS0lJgY2NTb55NjY2UKlUeP78eYHLhISEwMLCIm9ycHBQR1SdsOXqFjzLfobpPtOljlJmGto0xNYeW3FtzDW0cmqFLw59gQ8WfYDlF5dz8yaRjGlU2ZVGcHAwlEpl3pSQkCB1JNk4fu84mto3hV0FO6mjlDnnKs5Y33U94j6Ng7ejN0btG4W6v9TFmtg1yMnNkToeEZUxjSo7W1tbpKam5puXmpoKc3NzlCtX8DlTxsbGMDc3zzfR+xMEARH3IuBbw1fqKKJyruKMzd034/Koy3CzdcOQ3UPQaHkjHL97XOpoRFSGNKrsPD09ER4enm9eWFgYPD2lP4FZ19x4dAMPMh7Ivuxec7FxwY7eOxA9LBrlDcvj43Ufo9vWbrjz5I7U0YioDIhadunp6YiNjUVsbCyAV6cWxMbGIj4+HsCrTZBBQUF540eNGoW//voLkydPxo0bN/Drr79i27ZtGD9+vJgxqQAR9yJgoGcALwcvqaOolUdVD0R+EolN3Tbhwv0LqLe4HiaHTYYyUyl1NCJ6D6KW3YULF+Du7g53d3cAwIQJE+Du7o7p018d8JCcnJxXfADg5OSE/fv3IywsDK6urpg3bx5WrFgBf39/MWNSAS4lX4KLtYtO3lRVoVCgr0tf3Bh7A9NaTsPi84vzDmLh/jwi7cTLhVGB+mzvg0f/PMLRoKNSR5FckioJXx37Cusur0NDm4aY7z8fHzt9LHUsItmS/Xl2pDnSs9PL5N50clDVvCrWBq7N25/Xel1rdN/WHUmqJKmjEVExseyoQOnZ6Tq5CbMo/96fFxkfifq/1seS80t4zU0iLcCyowJlvMiAmSHL7k2v9+ddH3MdvT/sjU8PfIqWq1vi2sNrUkcjoiKw7KhA3IxZtIrlKmJ5p+WIGBSBh/88hNtSN8yMmMlLjxFpKJYdFYqb597Np4YPLo+6jCneU/D9qe/htswNp+NPSx2LiN7AsqMCVa1QFfef3Zc6hlYwMTDBrI9nIWZkDCxNLPHR6o8wat8opGWmSR2NiP4fy44KVM28GhJViVLH0CoNrBvg9JDTWNRuETZe2Yj6i+tjz809UsciIrDsqBAsu9LR19PHWI+xuPbpNTS2b4wuW7pg2J5heJb1TOpoRDqNZUcFqlqhKpLTk3nFkFJysHDAnj578Fun37Dl6ha4LXNDVELh92UkInGx7KhA1cyr4WXuSzzIeCB1FK2lUCgwrNEwxI6KhXV5a7RY3QLTjk3jffOIJMCyowJVM68GAIhXxr9jJL1L7Uq1cWrIKcz0mYmQ0yHwXOmJG49uSB2LSKew7KhAzlWcYahniIvJF6WOIgsGegaY5jMNUUOjkJ6djkbLGmHxucWQ2aVpiTQWy44KVM6wHBrZNUJkQqTUUWSladWmuDTyEoa4DcHYg2PRbmM7nuJBpAYsOyqUt4M3T5AWgamhKRZ3WIwD/Q7gcupluCxxwe4bu6WORSRrLDsqlLejN+KV8TwFQSTtPmiHK6Ov4CPHjxC4NRATD0/kwStEImHZUaG8HbwBAJHx3JQpliqmVbCz907M95+PhecWwmeNDxKUCVLHIpIdlh0VysbMBrUq1uJ+O5EpFAqMaz4Op4acQqIqEe7L3HHw1kGpYxHJCsuOiuTtyP126tK8WnPEjIxBs2rN0H5Te/wn/D94mftS6lhEssCyoyJ9XONjxKbE8ohBNalsWhl7++5FSOsQ/Bj5I/zW+SH5WbLUsYi0HsuOitS5bmcY6Bng92u/Sx1FZ+gp9DC1xVQcCzqGPx//Cbdlbjh295jUsYi0GsuOilSxXEX41fRD6LVQqaPoHJ8aPogZGQMXaxf4rfPDtye+5bVKiUqJZUfv1OvDXjgdf5qbMiVgY2aDwwMOY7rPdMyMmIlOmzvxPnlEpcCyo3fqUrcLN2VKSF9PHzN9Z+Jg/4OISoxC8xXNcfPRTaljEWkVlh29Ezdlagb/2v44N+wc9BR6aLaiGU9PICoBlh0VCzdlaoYPKn+As8POooVjC3TY1AE/Rf7Ei0kTFQPLjoqFmzI1h7mxOXb32Y2pLaZi8tHJCNoVhOcvnksdi0ijseyoWCqWq4h2H7TDypiVXJPQAPp6+vih9Q/Y3H0zfr/2O1quaYkkVZLUsYg0FsuOim1M0zG4nHqZlw/TIH0a9MHpT04jJT0FTX5rgrOJZ6WORKSRWHZUbH41/VCnch38cu4XqaPQvzSya4QLwy+gVsVa8FnjgzWxa6SORKRxWHZUbHoKPYxpOga/X/+dB6poGBszGxwbdAyDXAdhyO4hmHh4InKFXKljEWkMlh2VyCDXQTDWN8byi8uljkJvMNI3wrKOy7AwYCEWRC9Ar9BePHCF6P+x7KhELEwsEOQahGUXlyE7J1vqOPQGhUKBz5p9hp29d+LArQPwW++HR/88kjoWkeRYdlRiY5qOQUp6CnZc3yF1FCpE57qdETE4Arce34LXSi/ceXJH6khEkmLZUYl9aP0hWtVoxQNVNJxHVQ+cHXYWCoUCzVc2R3RitNSR3iIIwN69QHy81ElI7lh2VCpjPcYiMiESF+9flDoKFaFmxZo488kZ1K1cF63WtsKuG7ukjpTnzBmgZUtg8WLA0VHqNCR3LDsqlc51O6N2pdr4/tT3Ukehd6hsWhlHg46iY52O6La1GxZGL5Q0z7VrQGAg4O0NnD4N/PijpHFIR6il7BYvXowaNWrAxMQEzZo1w7lz5wodu2bNGigUinyTiYmJOmJSCRjoGSC4RTB23tiJqw+uSh2H3sHEwARbemzBRM+J+OLQF5hweILaT01ITASGDQNcXIDdu1/N698fcHNTawzSUaKX3datWzFhwgTMmDEDly5dgqurK/z9/fHgwYNClzE3N0dycnLe9Pfff4sdk0phYMOBcLRw5NqdltBT6OGntj9hUbtF+G/0f9V2asLTp8DUqcAHHwArVwK5/9+xhobArFmivz0RADWU3c8//4zhw4djyJAhqF+/PpYuXQpTU1OsWrWq0GUUCgVsbW3zJhsbm0LHZmVlQaVS5ZtIPQz1DTHVeyq2Xt3K+6tpkbEeY/NOTWi7oS2UmUpR3iczE5g7F6hVC5g9+9Xjfxs9GnByEuWtid4iatllZ2fj4sWL8PPz+98b6unBz88PUVFRhS6Xnp6O6tWrw8HBAV26dEFcXFyhY0NCQmBhYZE3OTg4lOnPQEX7xP0T2Fewx6yT/Ce6NulctzOODTqGuAdxaLW2FR5kFL6lpaRycoA1a4A6dYBJk16t2b2pQgXg66/L7C2J3knUsnv06BFycnLeWjOzsbFBSkpKgcvUrVsXq1atwu7du7Fhwwbk5ubCy8sLiYmJBY4PDg6GUqnMmxISEsr856DCGRsY46uPvsKmK5tw7eE1qeNQCTSv1hwnBp9AcnoyPlr9ERKUZfN3JysLMDEBmjR59d+CTJoEWFmVydsRFYvGHY3p6emJoKAguLm5wcfHBzt27ICVlRWWLVtW4HhjY2OYm5vnm0i9hroPhaOFI2ZGzJQ6CpWQi40LTg05heycbLRY3QJ/Pv7zvV/T1BTo0wdYvRpo2PDt521sgPHj3/ttiEpE1LKrUqUK9PX1kZqamm9+amoqbG1ti/UahoaGcHd3x+3bt8WISGXA2MAY01pOQ+i1UFxOuSx1HCqh2pVq4/SQ0yhvWB4frf4IsSmx7/2aSiXg7w8UdOD1jBmAmdl7vwVRiYhadkZGRmjcuDHCw8Pz5uXm5iI8PByenp7Feo2cnBxcuXIFdnZ2YsWkMhDkGoRaFWth2vFpUkehUqhqXhUnh5yEg7kDfNf4IjK+9PcsfF100f9/wRYrK2DBglf/X7v2q9MPiNRN9M2YEyZMwG+//Ya1a9fi+vXrGD16NDIyMjBkyBAAQFBQEIKDg/PGf/vttzhy5Aj++usvXLp0CQMGDMDff/+NYfwbotEM9Q0xq9Us7P1zL47dPSZ1HCqFKqZVcGzQMbjauqLthrY4cudIiV+joKI7dgwYMwaoUgX44YdXpxwQqZvoZde7d2/MnTsX06dPh5ubG2JjY3Ho0KG8g1bi4+ORnJycN/7p06cYPnw46tWrh/bt20OlUuHMmTOoX7++2FHpPfVp0Aee1Twx/vB45OTmSB2HSsHc2ByH+h9Cqxqt0HFTR2y/tr3YyxZWdA0aAAYGQEgI0KOHSMGJ3kEhCIIgdYiypFKpYGFhAaVSyYNVJHAu6RyarWiGZR2XYUTjEVLHoVJ6kfMCg3YNwta4rfit02/4xP2TIscXVXREJSXG97jGHY1J2s2jqgcGNhyIr499LdrJyiQ+Q31DrO+6HiMajcDQPUPx37P/LXQsi460AcuOytwPrX9AxosMXkZMy+nr6ePXDr9iktckjDs8DvPOzHtrDIuOtIWB1AFIfqqZV8MU7yn47uR3GNF4BGpXqi11JColhUKB2X6zYahniC/DvkSukItJ3pMAsOhIu3DNjkTxpdeXsDWzxaSwSVJHofekUCjw3cffYVrLaZh8dDJ+PP0ji460DtfsSBSmhqaY7Tcb/Xb0w7G7x/Cx08dSR6L3oFAo8G2rb6Gv0EdweDAWL8lBYvR/ALDoSDtwzY5E06dBHzSv1pynIsjIOPcZqPbnt0is8zXQ8jsWHWkNlh2JRqFQYIH/AvyR+gd+Pf+r1HHoPb3edJm4aRoQ/h3w8TT0+zWERUdagWVHompWrRlGNR6Fr459VWZX1Sf1e2sf3Y3/4NP63+C/cV9hTuQcacMRFQPLjkT3o9+PMDc2x6cHPoXMrmGgEwo7GGVxz+mY3nI6phydgrln5kobkugdeIAKic7CxAK/tPsF3bZ1Q+i1UPT6sJfUkaiY3nXU5UzfmcgRcjApbBL0FfoY78l795BmYtmRWnSt1xVdnbvi84Ofo03NNqhYrqLUkegdinN6gUKhwKxWs5CTm4MJRybASN8IYzzGSBOYqAjcjElq80v7X/D85XNMDpssdRR6h5KcR6dQKPBD6x8wvvl4jD04Fhv+2KDesETFwLIjtbGvYI/ZfrOxImYFTtw7IXUcKkRpThhXKBSY13YePnH7BIN3Dcaem3vUE5aomFh2pFYjGo9AC8cWGLFvBDJfZkodh97wPldGUSgUWN5pOQKdA9ErtBfva0gahWVHaqWn0MPyjstxL+0evjv5ndRx6F/K4hJg+nr62NhtI3xq+KDz5s44l3ROnLBEJcSyI7WrZ1UPX7X4CrMjZ+NK6hWp4xDK9qLOxgbG2NFrB1xtXdFuYztcfXC1bMMSlQLLjiQxtcVU1K1cF0G7gpCdky11HJ0mxkWdyxuVx/5+++Fg7oC269vir6d/lU1YolJi2ZEkjA2Msb7resQ9iMOM4zOkjqOzxLx7gaWJJQ4POAwzIzP4rfPD/Wf33/9FiUqJZUeScbdzx7etvsXsyNk49fcpqePoHHXcpsfGzAZHg47iZe5LtFnfBo//eVx2L05UAiw7ktQkr0nwdvTGwJ0DocpSSR1HZ6jzfnSOFo4IGxiGhxkP0W5jOzzLelb2b0L0Diw7kpS+nj7WBa7Dk+dP8MWhL6SOoxOkuPFq3Sp1cXjAYdx8fBNdt3blflpSO5YdSc6pohMWtluINbFrsOP6DqnjyJqUdxh3t3PHnj57cCr+FIbvHc6LgpNasexIIwxyHYSuzl0xYu8IJD9LljqOLElZdK/51PDBmi5rsO7yOsyMmKm+Nyadx7IjjfD66hsGegYYumco/9VfxjSh6F7r69IXP7b+Ed+e/BarYlapPwDpJJYdaYwqplWwqssqHLx9EEsuLJE6jmxoUtG9Ntl7MkY1HoURe0fg8O3D0gUhncGyI43S/oP2GN1kNL488iVuPLohdRytp4lFB7xak1/UfhECagegR2gPxKbEShuIZI9lRxpnbtu5cLRwRK/QXvjnxT9Sx9Famlp0rxnoGWBLjy2oW7kuOmzqgARlgtSRSMZYdqRxTA1NEdozFLef3MbnBz+XOo5W0vSie83MyAz7+u2DoZ4h2m9qD2WmUupIJFMsO9JILjYu+LXDr1gZsxJrY9dKHUeraEvRvWZrZouD/Q8iUZWI7tu68xw8EgXLjjTWYLfBGOI2BKP3j+aV84tJ24rutXpW9bC7z26eg0eiYdmRRvul/S+oVakWeob2RHp2utRxNJq2Ft1rLau3zDsHj/c6pLLGsiONZmpoiu09tyNRlYiR+0byX/yF0Paie62vS1984/sNpkdMx64bu6SOQzLCsiONV7dKXSzvuBybrmzCb5d+kzqOxpFL0b32dcuv0b1edwzcOZCbr6nMsOxIK/R16YvRTUbj84OfIyY5Ruo4GkNuRQcAego9rAlcg5oVa6Lz5s68LRCVCZYdaY2f/X/Gh9YfomdoTx6iDnkW3WtmRmbY3Wc3nmU/Q6/tvfAi54XUkUjLqaXsFi9ejBo1asDExATNmjXDuXPnihwfGhoKZ2dnmJiYwMXFBQcOHFBHTNJwJgYmCO0Zikf/PNL562fKueheq2FZA9t7bsfJv09i4pGJUschLSd62W3duhUTJkzAjBkzcOnSJbi6usLf3x8PHjwocPyZM2fQt29fDB06FDExMQgMDERgYCCuXuW2ewJqVqyJ1V1W4/frv+PH0z9KHUcSulB0r/nU8MHCgIVYdG4RVl5aKXUc0maCyDw8PIQxY8bkPc7JyRHs7e2FkJCQAsf36tVL6NChQ755zZo1E0aOHFng+MzMTEGpVOZNCQkJAgBBqVSW3Q9BGmfasWmCYqZC2HNjj9RR1O7TTwUBeDVZWQnClStSJxLfyL0jBcNvDYXTf5+WOgqpgVKpLPPvcVHX7LKzs3Hx4kX4+fnlzdPT04Ofnx+ioqIKXCYqKirfeADw9/cvdHxISAgsLCzyJgcHh7L7AUhjzfSdiS7OXdB/R39ce3hN6jhqFTLpCZqb/gErxUMcWxQnyzW6Ny1stxDNqzVHt23deA1NKhVRy+7Ro0fIycmBjY1Nvvk2NjZISUkpcJmUlJQSjQ8ODoZSqcybEhL4F0EX6Cn0sC5wHRwtHNFlSxc8ff5U6khqY/7TNBz+pwVOCh+hwagWwMWLUkcSnZG+Ebb32g4TAxMEbg3kBcKpxLT+aExjY2OYm5vnm0g3VDCugD199+DJ8yfovb03Xua+lDqSesyZA3OfRnDGTSAtDfDz04nCsy5vjd19duP6w+u8wACVmKhlV6VKFejr6yM1NTXf/NTUVNja2ha4jK2tbYnGk26rWbEmQnuG4tjdY5gcNlnqOOpRvjywfz/g4/PqsQ4VnputG1Z0XoENf2zAiksrpI5DWkTUsjMyMkLjxo0RHh6eNy83Nxfh4eHw9PQscBlPT8984wEgLCys0PFEHzt9jPn+8zH/7HzduUOCDhdeP5d+GNV4FD47+Blv+krFV2aHuhRiy5YtgrGxsbBmzRrh2rVrwogRIwRLS0shJSVFEARBGDhwoDB16tS88ZGRkYKBgYEwd+5c4fr168KMGTMEQ0ND4UoxDzkT4yge0ny5ubnC0N1DBeNZxsLZhLNSx1Gf9HRB8PH53+GZlpaCcOGC1KlE9/zFc8F9qbtQe2FtIe15mtRxqIyJ8T0ueriGc5YAACAASURBVNkJgiAsWrRIcHR0FIyMjAQPDw/h7Nn/fRn5+PgIgwYNyjd+27ZtQp06dQQjIyPhww8/FPbv31/s92LZ6a7MF5mC10ovwW6unZCkSpI6jvroaOHdfnxbMA8xF3ps6yHk5uZKHYfKkBjf4wpBkNdeXpVKBQsLCyiVSh6sooNS01PR5LcmsDOzw4nBJ1DOsJzUkdQjIwPo0AE4ceLVY0tL4OhRoHFjaXOJ7Pdrv6NHaA8sDFiIz5p9JnUcKiNifI9r/dGYRP9mY2aDXb134eqDqwjaFYRcIVfqSOqho/vwutfvji+afYGJRybiXFLRlyEk3cayI9lpbN8Ym7tvxu/XftedIzQBnS28OW3mwN3OHb1Ce+HJ8ydSxyENxbIjWeri3AUL2y3EvKh5WBS9SOo46qODhWekb4RtPbZBlaXC4F2Def4dFYhlR7I11mMsJnpOxBeHvsDuG7uljqM+Olh41S2rY13Xddj7517Mi5ondRzSQCw7krU5beage/3u6Pt7X0QnRksdR310sPA61umIKd5TMPXoVEQlFHwtXdJdLDuSNT2FHtZ3XY9Gdo3QaXMn3HlyR+pI6qODhffdx9+hadWmGLhzINKz06WOQxqEZUeyZ2Jggt19dsPSxBLtNrbDo38eSR1JfXSs8Az0DLC+63qkpKdg/KHxUschDcKyI51Q2bQyDvY/iLTMNHTZ0gXPXzyXOpL66Fjh1a5UG/P952NFzArsublH6jikIVh2pDNqVaqFvX33IiY5RrfOwQN0rvCGNRqGjnU6YtieYUhNT333AiR7LDvSKc2qNcs7B2/SkUlSx1EvHSo8hUKBFZ1e3RVh2N5hPB2BWHake16fg/fz2Z8x+/RsqeOolw4Vno2ZDVZ2Xol9f+7Db5d+kzoOSYxlRzpprMdYTG85HVPDp2LZhWVSx1EvHSq8TnU7YXij4Rh/eDxuPb4ldRySEMuOdNZM35n4zOMzjN4/GluubpE6jnrpUOH97P8z7MzsMHDnQN25mz29hWVHOkuhUGBBwAIMdB2IgTsHYv+f+6WOpF46UnhmRmZY33U9zt8/jx9O/SB1HJIIy450mp5CDys7r0THOh3RI7QHTv59UupI6qUjhefp4In/fPQffHviW1y8L6+fjYqHZUc6z0DPAJu7b4a3gzc6buqoe1+GOlJ401pOQwPrBhi+dzg3Z+oglh0RXl1lZVefXahvVR/+G/xx/eF1qSOplw4UnqG+IX7r9Bsup17G/Kj5UschNWPZEf0/MyMzHOh/AHYV7NBmfRvcS7sndST10oHCa1q1KT73+BwzImbgr6d/SR2H1IhlR/QvlcpVwpEBR2BiYAK/dX5ISU+ROpJ66UDhzfp4FqzLW2PUvlE82VyHsOyI3mBXwQ5Hg44i82Um2q5vq3t3v5Z54ZkZmeHXDr8i7K8wbPhjg9RxSE1YdkQFqGFZA2EDw5Ccnoy269vi6fOnUkdSL5kXXvsP2qNPgz4Yf3g8HmY8lDoOqQHLjqgQ9azq4ejAo7ibdhf+G/yRlpkmdST1knnhLfBfgFwhFxOOTJA6CqkBy46oCK62rggPCsftJ7cRsCEAqiyV1JHUS8aFZ2Nmg3lt52HDHxtw5M4RqeOQyFh2RO/gZuuGo0FHcfPxTQRsCMCzrGdSR1IvGRfeYLfB+NjpY4zaNwoZ2RlSxyERseyIiqGRXSOEDQzDtYfX0G5jOxaeTApPoVBgWcdlSE5Pxvenvpc6DomIZUdUTE3sm+DIwCO4+uAqAjZyk6ZcCq92pdqY5DUJP0f9jLtP70odh0TCsiMqAY+qHggbGIa4B3Hw3+APZaZS6kjqJdPCm+w9GZVNK2PK0SlSRyGRsOyISqhp1aY4GnQUNx7dQNsNbXmUpgwKz8zIDCGtQxB6LRSn/j4ldRwSAcuOqBSa2DdBeFA4bj2+hTbr2/A8PBkU3oCGA9DUvinGHR6HXCFX6jhUxlh2RKXUyK4Rjg06hrtP78JvvR8e/fNI6kjqJbPC01PoYUHAAlxKvoR1l9dJHYfKGMuO6D242brh2KBjSFQlouXqlkhSJUkdSb1kVnheDl7o06APgsODkZ6dLnUcKkMsO6L31NCmIU4NOYX07HS0WN0Cd57ckTqSesms8H5s/SPSMtMQcipE6ihUhlh2RGWgTuU6iPwkEsb6xmixugWupF6ROpJ6FafwFi0CcjV/X1h1y+r40vNLzIuap3u3eZIxlh1RGXGwcMDJISdha2YLnzU+iE6MljqSehVVeC9eADNmAJs3S5uxmKa0mIJK5SrxVAQZYdkRlSHr8tY4Pug4PrT+EK3XtUb4X+FSR1Kvwgrvp5+Ap0+Br78GsrKkzVgMZkZm+O7j77AtbhtikmOkjkNlQNSye/LkCfr37w9zc3NYWlpi6NChSE8veqevr68vFApFvmnUqFFixiQqU5Ymljg84DA+qv4R2m9qj103dkkdSb0KKrz//OfV/9+7ByxdKlm0kghyDcIHlT7A9IjpUkehMiBq2fXv3x9xcXEICwvDvn37cPLkSYwYMeKdyw0fPhzJycl505w5c8SMSVTmTA1NsbvPbnSp2wU9tvXA+svrpY6kXq8Lz9v77edmzQKUmn/lGQM9A8z0nYl9f+7TvU3SMiRa2V2/fh2HDh3CihUr0KxZM7Ro0QKLFi3Cli1bcP/+/SKXNTU1ha2tbd5kbm4uVkwi0RjpG2Fz980Y4jYEQbuCsCh6kdSR1CMjA1ixAujWDTh79u3nHz8G5s5Vf65S6P1hb9S3qs+1OxkQreyioqJgaWmJJk2a5M3z8/ODnp4eoqOL/lfSxo0bUaVKFTRo0ADBwcH4559/Ch2blZUFlUqVbyLSFPp6+ljeaTm+9PwSnx/6HLNOzIIgCFLHEpepKWBhAdy9C+TkFDzm55+B5GT15ioFfT19fOP7DY7cOYLT8aeljkPvQbSyS0lJgbW1db55BgYGqFSpElJSUgpdrl+/ftiwYQOOHz+O4OBgrF+/HgMGDCh0fEhICCwsLPImBweHMvsZiMqCQqHAnDZz8F2r7zA9YjpG7RuFl7kvpY4lHoUC6NkTiIsDliwBbGzeHvPPP8C336o/Wyl0q9cNDawb4LuT30kdhd6HUEJTpkwRABQ5Xb9+Xfj++++FOnXqvLW8lZWV8Ouvvxb7/cLDwwUAwu3btwt8PjMzU1AqlXlTQkKCAEBQKpUl/dGIRLfq0irB4FsDof3G9sKzrGdSx1GPZ88EYdYsQahQQRCA/036+oJw86bU6Yply5UtAmZCiE6MljqKTlAqlWX+Pa4QhJJtU3n48CEeP35c5JiaNWtiw4YNmDhxIp4+/d8Fcl++fAkTExOEhoaia9euxXq/jIwMmJmZ4dChQ/D393/neJVKBQsLCyiVSu7rI4105M4R9NjWAx9U/gD7+u6DXQU7qSOpx8OHwA8/AIsXvzrvDgB69ABCQ6XNVQw5uTn48NcPUadyHezpu0fqOLInyvd4mdXmG65duyYAEC5cuJA37/Dhw4JCoRCSkpKK/TqnT58WAAiXL18u1ngx/kVAVNZik2OFqvOqCtXnVxfiHsRJHUe9/vpLEAYMEASF4tUa3tmzUicqlrWxawXMhBCTHCN1FNkT43tctH129erVQ0BAAIYPH45z584hMjISY8eORZ8+fWBvbw8ASEpKgrOzM86dOwcAuHPnDmbNmoWLFy/i3r172LNnD4KCgtCyZUs0bNhQrKhEaudq64qzw87C3Ngc3qu8ceLeCakjqY+TE7B+PXDpEhAQAEyZ8mrDpobr59IPNSxr4Oeon6WOQqUg6nl2GzduhLOzM1q3bo327dujRYsWWL58ed7zL168wM2bN/OOtjQyMsLRo0fRtm1bODs7Y+LEiejevTv27t0rZkwiSVQzr4ZTQ06hiX0TtFnfBpuubJI6knq5uQEHD766jFhCgtRp3slAzwBjmo7B1riteJDxQOo4VEIl3men6bjPjrRNdk42RuwdgbWX1+KHj3/A1BZToVAopI5FBXjy/Amq/lwV01pOw1cffSV1HNkS43uc18YkkpiRvhFWd1mNGT4z8NWxr+R/aoIWq1SuEvq79MeSC0v4Z6RlWHZEGkChUGCm70ys6rwKq2JXofPmzlBl8QIJmmhM0zFIVCViz00elalNWHZEGmSI+xAc6HcAkQmRaL6iOW49viV1JHqDu507vBy8sPj8YqmjUAmw7Ig0TJtabRA9LBo5Qg48Vngg7E6Y1JHoDWObjsWxu8dw7eE1qaNQMbHsiDSQcxVnRA+LRvNqzRGwMQDzo+bL/5qaWqR7/e6wKW+Dxee4dqctWHZEGsrSxBL7+u7DRM+JmHBkAobsHoLMl5lSxyK8OqhoZOORWHt5LZ5lPZM6DhUDy45Ig+nr6WNOmzlY33U9tlzdAt81vkh+pvl3C9AFn7h/gowXGdj35z6po1AxsOyItMCAhgNwasgpJKgS0OS3JjifdF7qSDqvumV1eFT1QOg1zb+2J7HsiLRG06pNcWH4BThaOOKj1R9hwx8bpI6k83rW74kDtw5wU6YWYNkRaRG7CnaIGBSBfi79MHDnQEw6MoknN0uoZ/2eyMrJ4qZMLcCyI9IyxgbGWNl5JRb4L8D8s/Pht86P+/Ekwk2Z2oNlR6SFFAoFvmj+BY4POo5bT27BbZkbjt09JnUsndSrfi9uytQCLDsiLfZR9Y8QMzIGDW0awm+dH7498S1ycnOkjqVTetTvwU2ZWoBlR6TlrMtb41D/Q5jpOxMzI2ai3cZ2vAWNGr3elLn9+napo1ARWHZEMqCvp4/pPtMRNjAMl1Mvw32ZO079fUrqWDqjXe12iLgXgVwhV+ooVAiWHZGMtK7ZGjEjY1C7Um20WtsKcyLn8AtYDbwdvPHk+RPcfHRT6ihUCJYdkczYV7BHeFA4JntPxpSjU9BlSxc8ef5E6liy1rxac+gp9BCZECl1FCoEy45Ihgz0DPBD6x+wv99+nEk4w82aIqtgXAGuNq4sOw3GsiOSsfYftEfMyBg4mDvAZ40Pgo8GIzsnW+pYsuTt4I3IeJadpmLZEcmco4UjTgw+gR9a/4B5UfPQbEUzxD2IkzqW7Hg7euPWk1tITU+VOgoVgGVHpAP09fQxtcVURA+LRtbLLDRe3hgLzi7gwStlyNvBGwBwJuGMxEmoICw7Ih3ibueOiyMuYlSTURh/eDzarm+LRFWi1LFkwcHCAVamVrjy4IrUUagALDsiHVPOsBwWBCzAkQFHcP3RdbgsccGWq1ukjiULDhYOSFIlSR2DCsCyI9JRbWq1wZXRV+Bfyx99f++L/jv6Iy0zTepYWq2aeTUkPuOasiZi2RHpsErlKmFz983Y2G0j9v+5Hy5LXLD/z/1Sx9Ja1SpU45qdhmLZEek4hUKBfi798MfoP1Dfqj46bu6IPtv78KjCUqhqXpX7QDUUy46IALw6ReFQ/0NY33U9wu+Go97ielgVswqCIEgdTWtUM6+Gx88fI/NlptRR6A0sOyLKo1AoMKDhAFwfcx0d63TE0D1D0Xpda9x6fEvqaFqhaoWqAID7z+5LnITexLIjordUMa2CdV3X4fCAw7iXdg8uS1wQcioEL3JeSB1NKyigkDoCvYFlR0SFalurLa6MvoLPPD7D18e/RpPfmuBc0jmpY2ms9Ox0AICZkZnESehNLDsiKlJ5o/L4qe1POD/8PPQV+mi+ojk+P/g5nj5/KnU0jcOy01wsOyIqlkZ2jXBu+DnMaTMHq2NX44NFH2DxucV4mftS6mgaIz07HXoKPZgYmEgdhd7AsiOiYjPQM8CXXl/iz7F/okvdLvjs4GdwXeqKI3eOSB1NI6Rnp6O8YXkoFNxnp2lYdkRUYnYV7LCyy0qcH34elctVhv8Gf3Tc1FHn79T9+PljmBubSx2DCsCyI6JSa2zfGCcGn0Boz1DEPYxDgyUNMO7QOJ29M3pUYhSa2DeROgYVQLSy+/777+Hl5QVTU1NYWloWaxlBEDB9+nTY2dmhXLly8PPzw61bPL+HSJMpFAr0qN8D18dcx6xWs7AyZiU+WPQBfjn3i06dqpD5MhNRCVHwreErdRQqgGhll52djZ49e2L06NHFXmbOnDlYuHAhli5diujoaJQvXx7+/v7IzOTVCIg0nYmBCaa2mIo/x/6JwLqB+Pzg56i3uB7WX16PnNwcqeOJ7lzSOWTlZLHsNJRoZffNN99g/PjxcHFxKdZ4QRCwYMECfP311+jSpQsaNmyIdevW4f79+9i1a5dYMYmojL3enxczMgYNrBsgaFcQGixpgG1x22R9s9iIexGoaFIRDW0aSh2FCqAx++zu3r2LlJQU+Pn55c2zsLBAs2bNEBUVVehyWVlZUKlU+SYikp6rrSt29dmFc8POoYZlDfTe3htuS92wLW6b7Nb0BEHA7pu70bJ6S+gpNOZrlf5FY/5UUlJSAAA2Njb55tvY2OQ9V5CQkBBYWFjkTQ4ODqLmJKKSaVq1KQ72P4jITyJha2aL3tt7o8GSBtjwxwbZnKO34/oOXEq+hM+bfS51FCpEicpu6tSpUCgURU43btwQK2uBgoODoVQq86aEhAS1vj8RFY+XgxeODDyCqKFRqF2pNgbuHAjnX5yx5PySvCuPaKMXOS8QHB4M/1r++NjpY6njUCEMSjJ44sSJGDx4cJFjatasWaogtra2AIDU1FTY2dnlzU9NTYWbm1uhyxkbG8PY2LhU70lE6te8WnPs7bsXMckx+P7U9xh7cCymhk/FELch+LTpp6hTuY7UEUtkZcxK3H5yG6E9Q6WOQkUoUdlZWVnByspKlCBOTk6wtbVFeHh4XrmpVCpER0eX6IhOItIO7nbu2N5rO+KV8Vh6YSl+u/Qb/hv9X/jX8sdYj7FoV7sd9PX0pY5ZpARlAmZEzMCAhgPgausqdRwqgmj77OLj4xEbG4v4+Hjk5OQgNjYWsbGxSE//3+YKZ2dn7Ny5E8Crc3XGjRuH7777Dnv27MGVK1cQFBQEe3t7BAYGihWTiCTmaOGIH1r/gITxCVgbuBaPnz9Gp82d8MGiDzD3zFyNPUH9yfMnCNgYgHIG5TDbb7bUcehdBJEMGjRIAPDWdPz48bwxAITVq1fnPc7NzRWmTZsm2NjYCMbGxkLr1q2Fmzdvluh9lUqlAEBQKpVl9JMQkbpFJ0YLA3cMFIxmGQkm35kI3bZ2EzZf2Sw8y3omdTRBEAThvuq+4LrEVag8u7Jw4+ENqePIjhjf4wpBEATpqrbsqVQqWFhYQKlUwtyc16gj0mYPMh5g3eV12Ba3Defvn4eJgQnaf9AePev3RMc6HSW5lc7xu8cxdM9QZOdk49CAQ2hg3UDtGeROjO9xlh0RaYW7T+9i+7XtCL0Wmq/4etTrgVZOrWBrZivq+19OuYyp4VNx6PYhNKvaDNt6boOjhaOo76mrWHbFwLIjkr83iw8AalasCW8Hb7RwbAFvB2/Us6r33id4Z77MRHRiNFbGrMSGPzagdqXaCGkdgm71uvE2PiJi2RUDy45ItySpkhCZEInI+EhEJkQiNiUWOUIOKppUhKeDJ+pXqY+q5lVRzbwaqplXQ9UKVWFXwQ4Geq8ORhcEAc9fPkd6djrSs9MRr4zHiXsnEPF3BKISopCVkwU7MztM95mOoe5DYahvKPFPLH8su2Jg2RHptvTsdJxLOpdXfnee3kGSKgnPXz7PG6On0EOlcpWQ+TITGdkZEJD/a9DSxBI+1X3gW8MXrWq0gouNCy8DpkYsu2Jg2RHRmwRBwNPMp0hUJSJRlYgkVRIeZDxAOcNyMDMyg5mRGcobloeZkRmsy1ujvlV9jT/HT87E+B4v0UnlRETaSKFQoFK5SqhUrhLvSqCjuF5ORESyx7IjIiLZY9kREZHsseyIiEj2WHZERCR7LDsiIpI9lh0REckey46IiGSPZUdERLLHsiMiItlj2RERkeyx7IiISPZYdkREJHssOyIikj2WHRERyR7LjoiIZI9lR0REsseyIyIi2WPZERGR7LHsiIhI9lh2REQkeyw7IiKSPZYdERHJHsuOiIhkj2VHRESyx7IjIiLZY9kREZHsseyIiEj2WHZERCR7opXd999/Dy8vL5iamsLS0rJYywwePBgKhSLfFBAQIFZEIiLSEQZivXB2djZ69uwJT09PrFy5stjLBQQEYPXq1XmPjY2NxYhHREQ6RLSy++abbwAAa9asKdFyxsbGsLW1FSERERHpKo3bZxcREQFra2vUrVsXo0ePxuPHj4scn5WVBZVKlW8iIiL6N40qu4CAAKxbtw7h4eGYPXs2Tpw4gXbt2iEnJ6fQZUJCQmBhYZE3OTg4qDExERFpgxKV3dSpU986gOTN6caNG6UO06dPH3Tu3BkuLi4IDAzEvn37cP78eURERBS6THBwMJRKZd6UkJBQ6vcnIiJ5KtE+u4kTJ2Lw4MFFjqlZs+b75HnrtapUqYLbt2+jdevWBY4xNjbmQSxERFSkEpWdlZUVrKysxMrylsTERDx+/Bh2dnZqe08iIpIf0fbZxcfHIzY2FvHx8cjJyUFsbCxiY2ORnp6eN8bZ2Rk7d+4EAKSnp2PSpEk4e/Ys7t27h/DwcHTp0gW1a9eGv7+/WDGJiEgHiHbqwfTp07F27dq8x+7u7gCA48ePw9fXFwBw8+ZNKJVKAIC+vj7++OMPrF27FmlpabC3t0fbtm0xa9YsbqYkIqL3ohAEQZA6RFlSqVSwsLCAUqmEubm51HGIiKiExPge16hTD4iIiMTAsiMiItlj2RERkeyx7IiISPZYdkREJHssOyIikj2WHRERyR7LjoiIZI9lR0REsseyIyIi2WPZERGR7LHsiIhI9lh2REQkeyw7IiKSPZYdERHJHsuOiIhkj2VHRESyx7IjIiLZY9kREZHsseyIiEj2WHZERCR7LDsiIpI9lh0REckey46IiGSPZUdERLLHsiMiItlj2RERkeyx7IiISPZYdkREJHssOyIikj2WHRERyR7LjoiIZI9lR0REsseyIyIi2WPZERGR7LHsiIhI9kQru3v37mHo0KFwcnJCuXLlUKtWLcyYMQPZ2dlFLpeZmYkxY8agcuXKMDMzQ/fu3ZGamipWTCIi0gGild2NGzeQm5uLZcuWIS4uDvPnz8fSpUvx1VdfFbnc+PHjsXfvXoSGhuLEiRO4f/8+unXrJlZMIiLSAQpBEAR1vdlPP/2EJUuW4K+//irweaVSCSsrK2zatAk9evQA8Ko069Wrh6ioKDRv3vytZbKyspCVlZXvNRwdHZGQkABzc3NxfhAiIhKNSqWCg4MD0tLSYGFhUSavaVAmr1JMSqUSlSpVKvT5ixcv4sWLF/Dz88ub5+zsDEdHx0LLLiQkBN98881b8x0cHMomNBERSeLx48faV3a3b9/GokWLMHfu3ELHpKSkwMjICJaWlvnm29jYICUlpcBlgoODMWHChLzHaWlpqF69OuLj48vsQ1KX1/+a0ba1UuZWL+ZWP23Nrq25X2+hK2rlqKRKXHZTp07F7Nmzixxz/fp1ODs75z1OSkpCQEAAevbsieHDh5c8ZRGMjY1hbGz81nwLCwut+sP9N3Nzc63Mztzqxdzqp63ZtTW3nl7ZHVZS4rKbOHEiBg8eXOSYmjVr5v3//fv30apVK3h5eWH58uVFLmdra4vs7GykpaXlW7tLTU2Fra1tSaMSEREBKEXZWVlZwcrKqlhjk5KS0KpVKzRu3BirV69+Z0s3btwYhoaGCA8PR/fu3QEAN2/eRHx8PDw9PUsalYiICACgP3PmzJlivHBSUhJ8fX1RvXp1LFmyBM+fP0d6ejrS09NhZmaWN8bDwwMeHh6oWrUqTExMcP/+ffzyyy9wc3PDkydPMHLkSDg4OGDGjBnF/6H09eHr6wsDA7Uef1MmtDU7c6sXc6uftmZn7ldEO/VgzZo1GDJkSIHPvX7Le/fuwcnJCcePH4evry+AVyeVT5w4EZs3b0ZWVhb8/f3x66+/cjMmERGVmlrPsyMiIpICr41JRESyx7IjIiLZY9kREZHsseyIiEj2tL7stPlWQt9//z28vLxgamr61iXSCjN48GAoFIp8U0BAgMhJ8ytNbkEQMH36dNjZ2aFcuXLw8/PDrVu3RE76tidPnqB///4wNzeHpaUlhg4divT09CKX8fX1feszHzVqlKg5Fy9ejBo1asDExATNmjXDuXPnihwfGhoKZ2dnmJiYwMXFBQcOHBA1X2FKknvNmjVvfa4mJiZqTPvKyZMn0alTJ9jb20OhUGDXrl3vXCYiIgKNGjWCsbExateujTVr1ogf9A0lzR0REfHW561QKAq9FKNYQkJC0LRpU1SoUAHW1tYIDAzEzZs337nc+/6Oa33ZafOthLKzs9GzZ0+MHj26RMsFBAQgOTk5b9q8ebNICQtWmtxz5szBwoULsXTpUkRHR6N8+fLw9/dHZmamiEnf1r9/f8TFxSEsLAz79u3DyZMnMWLEiHcuN3z48Hyf+Zw5c0TLuHXrVkyYMAEzZszApUuX4OrqCn9/fzx48KDA8WfOnEHfvn0xdOhQxMTEIDAwEIGBgbh69apoGcsiN/DqMlb//lz//vtvNSZ+JSMjA66urli8eHGxxt+9excdOnRAq1atEBsbi3HjxmHYsGE4fPiwyEnzK2nu127evJnvM7e2thYpYcFOnDiBMWPG4OzZswgLC8OLFy/Qtm1bZGRkFLpMmfyOCzI0Z84cwcnJqdDn09LSBENDQyE0NDRv3vXr1wUAQlRUlDoi5rN69WrBwsKiWGMHDRokdOnSReRExVPc3Lm5uYKtra3w008/5c1LS0sTjI2Nhc2bN4sZMZ9r164JAITz58/nzTt48KCgUCiEpKSkQpfz8fERvvjiC3VEFARBEDw8PIQxY8bkPc7JyRHs7e2FkJCQAsf36tVL6NChQ755zZo1E0aOHClqzjeVNHdJfu/VBYCwIi8f7gAABvdJREFUc+fOIsdMnjxZ+PDDD/PN6927t+Dv7y9mtCIVJ/fx48cFAMLTp0/VlKp4Hjx4IAAQTpw4UeiYsvgd1/o1u4K8762ENF1ERASsra1Rt25djB49Go8fP5Y6UpHu3r2LlJSUfJ+3hYUFmjVrptbPOyoqCpaWlmjSpEnePD8/P+jp6SE6OrrIZTdu3IgqVaqgQYMGCA4Oxj///CNKxuzsbFy8eDHfZ6Wnpwc/P79CP6uoqKh84wHA399frZ9taXIDQHp6OqpXrw4HBwd06dIFcXFx6oj7XjTh834fbm5usLOzQ5s2bRAZGSl1HCiVSgAo8ju7LD5z7bp+TDGIdSshTREQEIBu3brByckJd+7cwVdffYV27dohKioK+vr6Uscr0OvP1MbGJt98dX/eKSkpb22yMTAwQKVKlYrM0a9fP1SvXh329vb4448/MGXKFNy8eRM7duwo84yPHj1CTk5OgZ/VjRs3ClwmJSVF8s+2NLnr1q2LVatWoWHDhlAqlZg7dy68vLwQFxeHatWqqSN2qRT2eatUKjx//hzlypWTKFnR7OzssHTpUjRp0gRZWVlYsWIFfH19ER0djUaNGkmSKTc3F+PGjYO3tzcaNGhQ6Liy+B3X2DW7qVOnFrgz9d/Tm3+JxLyVkJi5S6JPnz7o3LkzXFxcEBgYiH379uH8+fOIiIjQ6NxiEjv7iBEj4O/vDxcXF/Tv3x/r1q3Dzp07cefOnTL8KXSPp6cngoKC4ObmBh8fH+zYsQNWVlZYtmyZ1NFkqW7duhg5ciQaN24MLy8vrFq1Cl5eXpg/f75kmcaMGYOrV69iy5Ytor+Xxq7ZaeuthEqa+33VrFkTVapUwe3bt9G6detSv46YuV9/pqmpqbCzs8ubn5qaCjc3t1K95r8VN7utre1bB0u8fPkST548KdGfe7NmzQC82opQq1atEuctSpUqVaCvr//WkcFF/W7a2tqWaLwYSpP7TYaGhnB3d8ft27fFiFhmCvu8zc3NNXatrjAeHh44ffq0JO89duzYvIPE3rUmXxa/4xpbdtp6K6GS5C4LiYmJePz4cb4SKQ0xczs5OcHW1hbh4eF55aZSqRAdHV3iI1ELUtzsnp6eSEtLw8WLF9G4cWMAwLFjx5Cbm5tXYMURGxsLAO/9mRfEyMgIjRs3Rnh4OAIDAwG82tQTHh6OsWPHFriMp6cnwsPDMW7cuLx5YWFhar0tVmlyvyknJwdXrlxB+/btxYz63jw9Pd867F3dn3dZiY2NFeX3uCiCIOCzzz7Dzp07ERERAScnp3cuUya/46U9gkZTJCYmCrVr1xZat24tJCYmCsnJyXnTv8fUrVtXiI6Ozps3atQowdHRUTh27Jhw4cIFwdPTU/D09FRr9r///luIiYkRvvnmG8HMzEyIiYkRYmJihGfP/q+9uwdJLY7DOH6ofEGiTAi3IssaWqqhUYMgoiXachBpKKilhgSXiJoKooZotjEaIoegQKqhqKAXSEqiF6mppSEEbfve4dIhu/cOF67iPTwfOIvnJzwcfvyf5YhZc6atrY2trS0Astks09PTnJyckMlkSCaTdHV14ff7+fj4KNvcAAsLC7jdbhKJBNfX1wwODtLU1EQ+ny9ZboD+/n46Ozs5Ozvj6OgIv99PKBQy73/flYeHB+bn5zk/PyeTyZBIJPD5fAQCgaJl3NjYwOFwsL6+zu3tLWNjY7jdbl5fXwEIh8PEYjFz/vj4mKqqKpaWlkin08zOzmKz2UilUkXL+C9yz83Nsbe3x+PjIxcXFwwPD+N0Orm5uSlp7mw2a+6wYRgsLy9zdXXF8/MzALFYjHA4bM4/PT3hcrmIRqOk02nW1taorKxkd3e3rHOvrKywvb3N/f09qVSKyclJKioqSCaTJc09Pj5ObW0th4eHBed1LpczZ4qx4/992cXjcQzD+O31KZPJYBgGBwcH5mf5fJ6JiQnq6upwuVwMDQ0VFGQpRCKR3+b+mtMwDOLxOAC5XI6+vj7q6+ux2Ww0NjYyOjpqHiblmht+/vxgZmYGr9eLw+Ggt7eXu7u7kuYGeHt7IxQKUV1dTU1NDSMjIwUl/X1XXl5eCAQCeDweHA4HLS0tRKNR3t/fi5pzdXWVhoYG7HY73d3dnJ6emveCwSCRSKRgfnNzk9bWVux2O+3t7ezs7BQ135/8Te6pqSlz1uv1MjAwwOXlZckzf76S//36zBqJRAgGg798p6OjA7vdjs/nK9j1cs29uLhIc3MzTqcTj8dDT08P+/v7Jc/9p/P66zMsxo7rL35ERMTyyvZtTBERkX9FZSciIpanshMREctT2YmIiOWp7ERExPJUdiIiYnkqOxERsTyVnYiIWJ7KTkRELE9lJyIilqeyExERy/sB5p6g4dam49UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let us plot the ellipse along with the axes.\n",
    "# From our calculations, we know that the angle \n",
    "# of rotation is 45 degrees, and that the eigen vectors\n",
    "# are the columns of S\n",
    "\n",
    "import matplotlib\n",
    "fig = plt.figure(0)\n",
    "ax = fig.add_subplot(111, aspect='equal')\n",
    "\n",
    "e = matplotlib.patches.Ellipse((0, 0), 1, 3,\n",
    "                               theta, fc='None',\n",
    "                               edgecolor='g')\n",
    "# The ellipse is centered at (0, 0)\n",
    "# We are using random width and height. \n",
    "# Note that the direction of the axes is independent of\n",
    "# width and height\n",
    "ax.add_artist(e)\n",
    "ax.set_xlim(-2, 2)\n",
    "ax.set_ylim(-2, 2)\n",
    "ax.quiver([0, 0], [0, 0], S[:,0], S[:,1], color=['r','b','g'],\n",
    "          scale=5)\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
